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Foreword I 


Uncertainty is all pervasive—whether it relates to everyday personal choices and 
actions, or as background to business and policy decisions, or economic and climate 
predictions. In recent times, few things have attracted as much attention as the 
uncertainty surrounding the final whereabouts of MH370. 

How to deal scientifically with uncertainty? Put simply, on the one hand there 
are events or outcomes of interest that we don't know; on the other hand, pieces of 
information that we judge relevant in some sense that we do know. We need to 
assess what we believe about the unknowns, given the knowns. 

Formalising our measure of uncertainty in terms of probabilities, the scientific 
approach is encapsulated in the so-called Bayesian statistical paradigm, in which 
beliefs about the unknowns are quantified by a probability measure conditional on 
what we know. 

But typically, our state of knowledge itself gets modified over time and a method 
is therefore needed to refine and update beliefs as new information is acquired and 
assembled. The logical, mathematical rule for carrying out this updating is Bayes 
theorem, hence the term Bayesian Methods to describe the analytic and computa- 
tional toolkit that has been developed for updating beliefs as evidence changes or is 
added to. 

It is this toolkit that has been employed in the search for MH370 and this 
fascinating book provides a blow-by-blow case account of how the various strands 
of evidence have been brought together to give an overall probabilistic assessment 
of the final whereabouts of the plane. 

This has been an extremely complex task and the authors are to be congratulated 
on setting out systematically and coherently the science and mathematics driving 
the evidential equations. But, in addition to the complex modelling there remains 
the task of pulling out the Bayesian probability messages from the tangle of data 
that has been assembled. The computational methods for achieving this are of 
relatively recent origin and, on a personal note, I am delighted to have played a 
small part, with Dr. Neil Gordon, in the signal processing revolution that is now the 
particle filter method of analysis. 
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Dr. Gordon and his colleagues have employed these methods to provide a 
marvellous case study demonstrating the power of mathematical modelling and 
computation to attack one of the most intractable uncertain puzzles of recent times. 


London, UK Professor Sir Adrian Smith 
November 2015 FRS, Vice Chancellor 
University of London 


Foreword II 


The disappearance of the Malaysian Airlines flight MH370 from air traffic control 
radar on the evening of 8 March 2014 with 227 passengers and 12 crew on board 
was a tragedy and remains a mystery until we are able to locate the flight recorder or 
the wreckage. Only then can we unravel what actually occurred. 

Initially, the multinational search and recovery effort was focused on the Gulf of 
Thailand and the South China Sea. However, following closer analysis of Inmarsat 
data and flight path projections the search shifted to the Indian Ocean. 

The Malaysian government accepted the Australian government's offer to take 
the lead in the search and recovery operation in the southern Indian Ocean. 
Australia, led by the Australian Transport Safety Bureau (ATSB), took responsi- 
bility for defining the underwater search area in 28 April 2014. 

The Defence Science and Technology (DST) Group, being the only Australian 
government agency with the combined knowledge and experience in the scientific 
disciplines to support the search for MH370, became involved from May 2014. 
Defence scientists contributed a range of expertise across a spectrum of technolo- 
gies to the search: underwater acoustics, satellite communication systems and sta- 
tistical data processing. 

DST Group's position as a trusted government adviser with stewardship of the 
full range of defence technologies has been critical to our ability to contribute to the 
multi-agency—and multinational-search. Our ability to provide high quality, 
internationally respected, rapid response is built on our deep foundational research 
capability. Our support to the MH370 incident is an example of scientific research 
feeding directly into an active operational search. 

The quest for the MH370 demonstrates our ability to ‘work the full problem’ 
with a team of internationally-recognised experts contributing to the ATSB working 
group. Indeed, DST Group is proud to have been able to contribute its 
world-renowned expertise to the ATSB-led search. However, this book is not about 
the search for MH370. Rather, it focuses on the work to define the search zone. 

We are particularly fortunate to have Dr. Neil Gordon leading the DST Group 
team. Dr. Gordon is recognised internationally as an expert in statistical data 
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processing and, in particular, the dynamic Bayesian estimation methods deployed in 
this activity. 

Working alongside Dr. Gordon are other Defence scientists—Dr. Samuel Davey, 
Dr. Mark Rutten and Dr. Jason Williams who are all experts in target tracking and 
multi-sensor fusion, and Dr. Ian Holland who specialises in satellite and wireless 
communications. 

DST Group's expertise and unique capabilities will continue to contribute to the 
ongoing search for the MH370, led by the ATSB and in collaboration with other 
Australian and international agencies. 


Canberra, Australia Dr. Alex Zelinsky 
November 2015 Chief Defence Scientist 
Defence Science and Technology 

Department of Defence 
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Chapter 1 
Introduction 


On 7 March 2014 at 16:42,! Malaysian Airlines flight MH370 departed from 
Kuala Lumpur (KL) International Airport bound for Beijing. There was a total 
of 239 persons on board (227 passengers and 12 crew). The aircraft was a 
Boeing 777-200ER registered as 9M-MRO. The aircraft lost contact with Air 
Traffic Control (ATC) during a transition between Malaysian and Vietnamese 
airspace. The last recorded radio transmission from MH370 was at 17:19. Over 
the following days, an intensive air and sea rescue search was made around the 
last reported position of the aircraft in the Gulf of Thailand without success. It 
then became clear that satellite communication messages between the aircraft 
and one member of Inmarsat's constellation of geosynchronous satellites were 
crucial to defining the search zone for the aircraft. Satellite communication systems 
involve transmissions between multiple terminals using a satellite. The aircraft 
9M-MRO communicated with the Inmarsat ground station in Perth, Western 
Australia, via the Indian Ocean Region satellite Inmarsat-3F1. The data available for 
MH370 is mostly comprised of approximately hourly “handshake” transmissions 
initiated by the ground station for aircraft that have not communicated in the 
preceding hour. No explicit information relating to the aircraft terminal location is 
contained in the messages; however, the messages contain metadata which can be 
processed to produce estimates of the flight path and final location. 

Inmarsat conducted a rapid and innovative analysis of the data which placed the 
aircraft in the Australian search and rescue zone on an arc in the Southern part of the 
Indian Ocean. On 17 March 2014, the Australian Maritime Safety Authority took 
responsibility of the search and rescue operation. Subsequently the Joint Agency 
Coordination Centre (JACC) was established on 30 March 2014 to coordinate the 
Australian Government's support for the search and over the following weeks an 
intensive aerial and surface search was conducted by an international team. 


‘All times are given in Coordinated Universal Time (UTC) in format hh:mm:ss. Local time in 
Malaysia and Western Australia is 8h ahead of UTC and local time in West Indonesia is 7h ahead 
of UTC. 
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On 28 April 2014, the aerial search concluded and the search moved to an 
underwater phase. The Australian Transport Safety Bureau (ATSB) took respon- 
sibility for defining the underwater search area. The ATSB convened a flight path 
prediction working group in order to bring together experts in satellite communi- 
cation systems and statistical data processing and apply novel data analysis tech- 
niques to estimate the most likely final location of MH370. This working group 
consisted of representatives from the following organisations: Air Accidents Inves- 
tigation Branch (UK); Boeing (US); Inmarsat (UK); National Transportation Safety 
Board (US); Thales (UK); and the Defence Science and Technology (DST) Group 
(Australia). The working group developed new methods to analyse the Inmarsat data 
and validated those methods. The ATSB released a report summarising the findings 
of the working group in August 2014 [3] and have subsequently released updates in 
October 2014 [4] and December 2015 [5]. 

In this book we detail the statistical approach adopted by the DST Group team 
to analyse the available data and produce a probability density function (pdf) of 
the accident aircraft's final location. In Chap. 2, we start by detailing a summary of 
the events and listing of the available data. In Chap.3, we describe the Bayesian 
framework, upon which our method is formulated. This approach requires several 
ingredients. The first of these is a prior distribution; this is built on the primary 
radar data described in Chap. 4. The second ingredient consists of a set of likelihood 
functions detailing how the available measurements are linked to the aircraft state; 
the measurement models and error statistics are characterised in Chap. 5. The final 
ingredient is a stochastic model describing the possible dynamic trajectories of the 
aircraft. Chapters 6 and 7 describe our models for aircraft dynamics during cruise and 
manoeuvre respectively. The resulting set of (prior, likelihood, dynamics) enables 
us to calculate the probability distribution of aircraft trajectories. However since the 
models are nonlinear and non-Gaussian we are required to use numerical methods for 
the calculation; our particle filtering approach is described in Chap. 8. The Bayesian 
method that was developed has been validated against a number of earlier flights of 
the accident aircraft 9M-MRO, where accurate measurements of the aircraft loca- 
tion were available from the aircraft's logging system. These results are detailed in 
Chap. 9. The method has also been applied to the data available for the accident flight; 
the resulting probability distribution, which defines the search zone, is described in 
Chap. 10. The search zone is defined by combining our pdf from the analysis of 
the satellite data with a kernel describing the distribution of aircraft motion during 
descent, which was defined by expert accident investigators from the ATSB. Thus, 
any adjustment to the assumptions made about the descent (and hence the kernel 
describing its distribution) yields a change to the search region. Finally, in Chap. 11, 
we discuss on-going work including the impact of the flaperon wreckage discovery 
on Reunion Island. 

Being statistical in nature, this book seeks to identify aircraft paths that are most 
likely, given models which incorporate sensor measurements, commonly observed 
commercial aircraft motion, and air transport safety investigator and manufacturer 
expert advice and assessments. The goal is not to identify the complete and exhaus- 
tive set of all possible aircraft paths, but rather the subset of those trajectories that 
are most probable. The area covered by the set of end-points of all possible paths is 
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prohibitively large and the overwhelming majority of this area contains vanishingly 
small probability. In determining how to allocate finite search resources, priority 
is directed to the area containing the highest probability trajectories. The proposed 
method incorporates a rich variety of possible paths, and calculates a probability 
distribution based on the well-established, rigorous Bayesian toolkit, which auto- 
matically trades the complexity of the model against the match to the observed 
measurements. 

The ATSB maintains a website http://www.atsb.gov.au/mh370.aspx with compre- 
hensive information about the ongoing search and we refer readers to this resource 
for the latest news, data, maps, videos, reports and operational updates. Any com- 
ments or information related to the search or the analysis in this book should be sent 
to atsbinfo @atsb.gov.au. Feedback related to the analysis in this book can also be 
sent to MH370QG dsto.defence.gov.au. 


1.4 Summary of Results 


The key outcome of this book is a probability distribution of the final aircraft 
location based on a Bayesian analysis, using models constructed through detailed 
study of the measurement noise statistics and commercial aircraft motion, incorpo- 
rating assessments of likely operating parameters from expert accident investigators. 
The distribution is shown in Fig. 1.1, illustrated as a contour plot, where the colour 
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Fig. 1.1 Probability density function of final location of MH370. Indicative search area (as of 
November 2015) marked with solid black line 
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Google earth 
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Fig. 1.2 Result of validation analysis applied to the 9M-MRO flight from KL to Amsterdam on 26 
February 2014 


indicates the likelihood, with red being most likely, and blue being least likely. The 
black rectangle shows the indicative search area as at November 2015. Full details 
can be found in Chap. 10. In order to orient the scale of the search being undertaken, 
the area of the indicative search region shown in Fig. 1.1 is 100,000 km?. In contrast, 
the bounding region considered in the search for Air France flight 447? was a disc 
with radius 40 nm [40], corresponding to an area of 17,000 km?. 

The procedure used to generate the result in Fig. 1.1 was validated by applying the 
identical process to other flights for which the location was known. The flights exam- 
ined included four flights for the same aircraft (OM-MRO) involved in the accident, as 
well as two others that were in flight at the same time as the accident. The communi- 
cations logs for these flights were down-sampled to emulate the information available 
on the accident flight. An example of the result is shown in Fig. 1.2. The white line 
shows the true path taken by the aircraft, while other coloured lines show potential 
trajectories, coloured by their probability, again with red indicating the most likely 
paths and blue indicating the least likely. The probability distribution at the time of 
the final measurement is shown as a red line, with the likelihood encoded through 
the height above the earth. 

Ten different subsets of measurements were used from each of the six flights 
to make a total of 60 validation experiments. In all cases, the true final aircraft 


?Air France Flight 447 disappeared over the South Atlantic Ocean on 1 June 2009. The initial 
underwater search consisted of three phases over approximately one year. After the third phase of 
the search, US company Metron conducted a Bayesian analysis of the available data. Debris was 
found in April 2011 close to the final reported location of the aircraft. 
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location was contained within the 85 ?6 probability region. This indicates that the 
probability distribution produced by the method is slightly conservative, otherwise 
only around 50 of the 60 experiments should be in the 85 % region. This is because 
the dynamics model allows for a wider range of aircraft manoeuvres than are 
actually experienced by typical commercial flights. Such a model is appropriate 
because the accident flight was not a typical commercial flight. Full details of the 
validation experiments can be found in Chap. 9. 


Open Access This chapter is distributed under the terms of the Creative Commons 
Attribution-NonCommercial 4.0 International License (http://creativecommons.org/licenses/ 
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Chapter 2 
Factual Description of Accident 
and Available Information 


The detailed chronological factual statement of known information about flight 
MH370 is given at [35]. A brief summary is given here sufficient to put the analysis 
in the rest of the book in context. 

On 7 March 2014 at 16:42, flight MH370 departed from KL bound for Beijing; 
this is marked as event 1 in Fig. 2.1. Initially, everything about the flight proceeded 
as normal. The Mode S transponder system on-board the aircraft was responding as 
expected to interrogation from the ATC Secondary Surveillance Radar up to the time 
when it was lost on the ATC radar screen at 17:21:13, marked as event 4. No mes- 
sage was received from the aircraft to report a system failure. Similarly the on-board 
Aircraft Communications Addressing and Reporting System (ACARS) reported as 
expected at 17:07:29, this event is marked as 2 in Fig. 2.1. This message contained a 
collation of six reports generated at five minute intervals by the system from 16:41:43 
until 17:06:43. These reports contained information about the aircraft position and 
motion such as latitude, longitude, altitude, air temperature, air speed, wind direc- 
tion, wind speed, and true heading. The ACARS position reports are scheduled to 
be transmitted at thirty minute intervals during cruise. The next scheduled report at 
17:37 was not received. The last recorded radio transmission with the crew of MH370 
occurred at 17:19:30 as the aircraft was instructed to contact Vietnamese ATC on 
leaving the Malaysian Flight Information Region. This communication is marked 
as event 3. At 17:39:06 Vietnamese ATC contacted KL ATC to query the where- 
abouts of MH370, who initiated the efforts of several countries' ATC to establish the 
location, without success. Malaysian military radars were subsequently able to show 
primary radar returns associated with MH370 deviating from the flight plan almost 
immediately after the loss of Secondary Surveillance Radar at 17:21:13 by making 
a left turn to end up travelling in a South Westerly direction. Radar returns show 
the aircraft travelling back across Malaysia before turning near Penang Island and 
travelling in a North Westerly direction up the Straits of Malacca. The final primary 
radar return was recorded at 18:22:12, which is marked as event 5 in Fig.2.1. 
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‘Kuala Terengganu 


Fig. 2.1 Chronological sequence of events in MH370 disappearance. / = Take off; 2 = Final 
ACARS transmission; 3 = Final radio contact; 4 = Final Mode S transmission; 5 = Final primary 
radar detection 


The ACARS system also provides message communication between the aircraft 
and its ground base. At 18:03:23 Malaysian Airlines Operations Centre sent a mes- 
sage asking the crew to contact Vietnamese ATC immediately. The downlink message 
showed that this failed to reach the aircraft. The message was auto transmitted every 
two minutes until 18:43:33 with no success [35]. The next communication on the 
SATCOM system was a log-on request from the aircraft at 18:25. This communica- 
tion contains no location information from the aircraft but it does however contain 
timing and frequency metadata which turn out to be crucial in estimating where 
MH370 went after the final radar return at 18:22:12. The timing measurement is 
called the Burst Timing Offset (BTO) and the frequency measurement the Burst 
Frequency Offset (BFO). At 18:39:52 there was an unanswered attempt to call the 
aircraft from ATC on the ground to air telephone link. The failed communication 
data log contains frequency but not timing metadata. 

If there has been no SATCOM activity for sixty minutes then the ground station 
automatically initiates a handshake confirming presence of the aircraft. If the air- 
craft receives this query then an automatic response is sent indicating that it is still 
logged on to the SATCOM network. These handshake responses contain both tim- 
ing and frequency metadata. Handshake messages occurred at 19:41:00, 20:41:02, 
21:41:24 and 22:41:19. At 23:15:58 there was a second unanswered attempt at 
a satellite telephone call, again giving frequency metadata. There was a ground 
station initiated handshake at 00:10:58 (timing and frequency metadata). Finally, 
at 00:19:29 the aircraft SATCOM system initiated another log-on request and this 
was the last SATCOM transmission received from the aircraft. Due to its timing, 
this log-on request is believed to correspond to fuel exhaustion and subsequent 
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Table 2.1 Summary of SATCOM data available for MH370 


Event Time (UTC) BTO BFO 
Aircraft departed KL 16:42 Y Y 
Last ACARS transmission 17:07 Y Y 
AES initiated log on 18:25 Y N* 
AES access request 18:28 Y Y 
Unanswered ground to air telephone call | 18:39 N Y 
GES initiated handshake 19:41 Y Y 
GES initiated handshake 20:41 Y Y 
GES initiated handshake 21:41 Y Y 
GES initiated handshake 22:41 Y Y 
Unanswered ground to air telephone call | 23:15 N Y 
GES initiated handshake 00:10 Y Y 
AES initiated logon 00:19 Y N* 


Measurements marked with an asterisk are available but cannot be used as discussed in the text 


activation of the auxiliary power unit [5]. There was no response to the ground 
station initiated handshake request at 01:15:56. The log-on messages at 18:25 and 
00:19 contain BFO measurements but the equipment was not in steady state. The 
values of the BFO measurements were deemed to be unreliable at these times and 
cannot be used. 

The key pieces of information available to us to estimate the MH370 flight path 
are the final radar detection at 18:22:12 and the timing and frequency metadata 
associated with the infrequent SATCOM messaging that subsequently occurred up 
until 00:19. Table 2.1 summarises the available measurement data. The messages at 
18:25 and 00:19 occurred during transient phases of operation for the SATCOM 
equipment so the BFO values reported for these times cannot be used. In the 
next two chapters we describe the Bayesian approach, and then how the satellite 
communication system works and how the timing and frequency metadata can 
be related to aircraft location, allowing us to build a likelihood function for our 
Bayesian analysis. 
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Chapter 3 
The Bayesian Approach 


Bayesian inference methods [9] provide a well-studied toolkit for calculating 
a distribution of a quantity of interest given observed evidence (measurements). 
As such, they are well-suited for calculating a probability distribution of the final 
location of the aircraft given the data available from the Inmarsat satellite 
communication system. The resulting probability distribution is essential to 
prioritise search efforts. In this chapter, we provide a brief introduction to Bayesian 
methods. We assume a reasonable background in probability theory; the interested 
reader is referred to excellent resources such as [8, 9, 10, 18, 36, 39, 43] if further 
detail is desired. 

The required probability density function (pdf) is the probability of the aircraft 
location given the available data. Bayes' rule defines a method to calculate this pdf 
using prior information, including knowledge of how aircraft move, and a model of 
how the measured data relate to the aircraft location and velocity. Mathematically, 
Bayes' rule is 


pxiz) = ae G.D 
i — (3.2) 
(0 peal) p) = 


J pix?) p(x’)dx’ 
where the elements are: 


1. x is the random variable, or the state, which is the quantity of interest (e.g., the 
position of the aircraft); 

2. z is the measurement (e.g., the Inmarsat satellite communication data, which 
provides some form of positional data); 

3. p(x) is the prior pdf of the state (not incorporating the measurement, e.g., based 
on historical data); 
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4. p(z|x) is the pdf of the measurement conditioned on the state (e.g., this may be 
constructed by observing the distribution of measurements in cases where the 
state is known); 

5. p(x|z) is the conditional pdf of interest (the posterior pdf), describing the distribu- 
tion of state (e.g., aircraft location) taking into account the observed measurement. 


The posterior probability density is based on the accumulated Inmarsat satellite 
communications data as well as all available contextual knowledge on the sensor 
characteristics, aircraft dynamic behaviour and environmental conditions and con- 
straints. The method is based on the state space approach to time series modelling. 
Here, attention is focused on the state vector of a system. The state vector contains 
all relevant information required to describe the system under investigation at a given 
point in time. For example, in radar tracking problems this information would typi- 
cally be related to the kinematic characteristics of the aircraft, such as position, 
altitude, speed, and heading. The measurement vector represents noisy observations 
that are related to the state vector. For example, the distance and bearing angle 
between the sensor and the object being measured. The state-space approach is 
convenient for handling multivariate data and nonlinear, non-Gaussian processes; 
it provides a significant advantage over traditional time series techniques for these 
problems; and has been extensively used in many diverse applications over the last 
50 years [7]. An excellent summary of Bayesian techniques for state space models 
is given by [36]. 

In order to proceed, two models are required: first, the measurement model relates 
the noisy measurements to the state; and second, the system or dynamic model 
describes the evolution of the state with time. The measurement model used for BTO 
and BFO metadata is defined in a probabilistic form in Chap. 5. The dynamic model 
used to define the behaviour of the aircraft is defined in Chaps. 6 and 7. 

If the measurement model and the system model are both linear and Gaussian, 
the optimal estimate can be calculated in closed form using the Kalman filter [25]. 
If either the system or measurement model is nonlinear or non-Gaussian, the pos- 
terior pdf will be non-Gaussian and standard analysis with a Kalman filter will 
be suboptimal. This results in the need for approximate computational strategies 
and the approach adopted in this study is introduced in this chapter. The appli- 
cation of the measurement and dynamics models to this approach is described in 
Chap. 8. The computational approach proceeds in essentially two stages: prediction 
and update. The prediction stage uses the aircraft dynamic model to step from the 
state pdf at one time to the pdf at the next time. The state is subject to unknown 
disturbances, modeled as random noise, and also unknown control inputs, such as 
turn commands, and so prediction generally translates, deforms, and broadens the 
state pdf. The update operation uses the latest measurement to modify (typically to 
tighten) the prediction pdf. This is achieved using Bayes theorem, (3.3), which is the 
mechanism for updating knowledge about the state in the light of extra information 
from new data. 
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3.1 The Problem and its Conceptual Solution 


To define the problem of nonlinear filtering, let us introduce the state vector 
x(t) € R”, where n is the dimension of the state vector. Here ¢ is continuous-valued 
time. The state evolution is best described using a continuous-time stochastic differ- 
ential equation, sometimes specifically referred to as an Itó differential equation [23]. 
However, itis often more convenient to sample this at discrete time instants, in which 
case Xy = X (tg) represents the state at the kth discrete sample time. The elapsed time 
between samples A, = tk — tk—1 is not necessarily constant. The state is assumed to 
evolve according to a continuous-time stochastic model: 


dx(t) = f (x(t), dv(t), t, dt), (3.4) 


where f(-) is a known, possibly nonlinear deterministic function of the state and v(t) 
is referred to as a process noise sequence, which caters for random disturbances in 
the aircraft motion. 

A sensor collects measurements, which are a possibly nonlinear function of the 
state. Measurements occur at times tg, for k = {1,2,...K}. The kth measurement 
is denoted z; € IR" where m is the dimension of the measurement vector. The mea- 
surements are related to the state via the measurement equation: 


Zk = h; (Xk, We) (3:5) 


where h, (-) is a known, possibly nonlinear function and w; is a measurement noise 
sequence. The noise sequences v(t) and w; will be assumed to be white, with known 
probability density functions and mutually independent. The initial state is assumed 
to have a known pdf p (xo) and also to be independent of noise sequences. 

We seek estimates of x, based on the sequence of all available measurements up to 
time f;, defining the measurement history Z, = (zi, .. . Z}. From a Bayesian perspec- 
tive, the problem is to recursively construct the posterior pdf p (x;|Z;). In principle, 
the pdf p (x;|Z;) may be obtained recursively in two stages: prediction and update. 
The prediction stage steps from the pdf of x at time fj 1, p (Xi 1|Z41), to the pdf 
at the next time, p (x;,|Z,_;), not incorporating any new measurements. The update 
stage takes the predicted pdf p (x,|Z_1) and incorporates the new measurement Zg 
occurring at time f to obtain the updated pdf p (x;,|Z;). If there is a requirement 
to evaluate the pdf at time ¢ for which there is no measurement then this pdf is the 
predicted pdf and no update step needs to be performed. 


3.1.1 Prediction 


The prediction stage involves using the system model (3.4) to obtain the prediction 
density of the state at time step k via the Chapman-Kolmogorov equation: 
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pez.) = | »eux 1, Zia) p Ga-alZic) AXK-1, 


ES 1) P Gu 1Zia) Ax. (3.6) 


The first line of (3.6) is a statement of the law of total probability. The sim- 
plification p (X.|Xx 1, Zi 1) = p (Xk|Xz—1) used to progress from the first line of 
(3.6) to the second applies because (3.4) describes a Markov process of order one. 
The probabilistic model of the state evolution, p (x; |Xi.1), is defined by the system 
equation (3.4) and the known statistics of v(t). 


3.1.2 Update 


At time f, a measurement z; becomes available and the update stage is carried out. 
This involves an update of the prediction (or prior) pdf via Bayes' rule: 


P 074) = p Quz, Zia) 
P (2x|Xx, Zi) p Gud) 
P GdZica) 
Z- 
= P G|Ixi) p Gu|Z D (3.7) 
P (2x |Ze_1) 


where conditional independence has been used to write the likelihood function 
P (2x |X, Z1) = p (Z |X), which is defined by the measurement model (3.5) and 
the known statistics of w;. The normalizing constant on the denominator can be 
expanded as 


P (Zk|Zk-1) = J pais) p (Xk |Z—1) AXx. (3.8) 


In the update stage (3.7), the measurement z; is used to modify the prior density to 
obtain the required posterior density of the current state. 

Note that there is no requirement for all of the measurements to have the same 
statistical model or even contain the same type of information. For example, there 
could be multiple sensors operating on different modalities. For simplicity, we have 
not introduced explicit notation to change the measurement pdf for each k. For the 
accident flight three different types of measurement have been used. As discussed in 
Chap. 5, the satellite communications messages consist of R-channel and C-channel 
messages that have differing information content. Another quite different form of 
measurement is the areas of the ocean floor that have been searched without locating 
the aircraft and the debris that has been recovered. This measurement and its potential 
use to refine the ongoing search are discussed in Chap. 11. 
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The recurrence relations (3.6) and (3.7) form the basis for the optimal Bayesian 
solution. The recursive propagation of the posterior density, given by (3.6) and (3.7), 
is only a conceptual solution in the sense that in general it cannot be determined 
analytically. In most practical situations the analytic solution of (3.7) and (3.8) is 
intractable and numerical approximations have to be used. This has been a topic of 
significant research effort over the past 20 years [1, 20, 33]; a general overview of 
the method is presented next. 


3.20 The Particle Filter 


In the linear Gaussian case, the pdfs for p (v(t)), p (Wx), p (Xo) are all Gaussian and 
the functions f(-) and h(-) are linear. It can then be easily shown that the posterior 
P G| Z4.) is also Gaussian and all of these pdfs can be summarised completely by 
their means and covariances. The Kalman filter is an algorithm that defines recursions 
for the mean and covariance of p (x;|Z,) in terms of the means and covariances of 
the prior and noise processes. However, in general, the posterior does not take the 
same functional form as the prior and indeed it is not possible to even write a closed 
form expression for p (x;|Z;). In this case an approximate solution is required. 
The solution used for the MH370 search definition is referred to as the particle filter 
and is a numerical approximation based on random sampling. 

The fundamental concept in the particle filter is to approximate the pdf p (x;|Zx) 
as a weighted combination of sample points 


P 
p lZ ~ S wed (x, — x). (3.9) 


p=1 


where the w? are referred to as weights and sum to unity, and the x? are referred 
to as particles. The convergence properties of this approximation in the limit as the 
number of particles P increases have been well studied, for example [14, 21]. Given 
this approximate pdf, itis simple to evaluate the expectation of any nonlinear function 
of the state, such as 


P 


[eo iz] = [aeo p euZocx Y wed). G0 


p=1 


The approximation of an integral using sample points as above is referred to as 
Monte Carlo integration and can be applied to both the Chapman—Kolmogorov pre- 
diction (3.6) and the Bayesian update (3.7). 

The particle filter is an algorithm that provides a mechanism to recursively create 
a set of weighted particles approximating p (x;|Z;,) starting from a previous set of 
weighted particles approximating p (Xi 1|Z4—1). It does this in two stages: first it 
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moves the particle sample points x? , — x; to new locations using a pdf referred to 
as a proposal distribution, which is a tractable approximation of the pdf of interest. 
Second, it determines new particle weights to correct for the difference between the 
proposal and the true pdf. This process is known as importance sampling [1, 33]. 

The proposal distribution is a critical component of the particle filter. It is a 
function chosen by the designer subject to relatively loose constraints. Importantly, 
the proposal distribution must cover all of the state space where the true distribution 
is non-zero and its tails should be heavier than the tails of the true distribution. If the 
proposal is chosen poorly then many of the particles x will be assigned very low 
weights and the filter efficiency will be low: a large number of particles will be 
required for satisfactory performance. A common version of the particle filter is the 
Sample-Importance-Resample (SIR) particle filter that uses the system dynamics as a 
proposal distribution. The SIR is popular because it is often relatively straightforward 
to sample from the dynamics and because the weight update equation is very simple 
when the dynamics is used as the proposal. The filter used in this book is a form of 
SIR particle filter. 

For the SIR particle filter, for each particle x , a new x; is drawn from the tran- 
sition density p(x,|x7. ;), and weights are updated by scaling the previous weights 
by the current measurement likelihood and re-normalising, 


we = (Wi! p (zidxt) wr 4. (3.11) 


where the normalising term is Wy = X p (dx) wL 4. 

A key difficulty in particle filters is the issue of degeneracy, i.e., over time, many 
weights tend toward zero, and the corresponding particles are of little use. Resam- 
pling is used to combat this difficulty. The simplest approach is draw P new particles 
from the approximate distribution (3.9), such that particles with very large weights 
are likely to be replicated many times over, and those with very small weights 
are unlikely to be sampled. A variety of methods are possible, and can be found 
in [1, 33]. The sampling method used in this study is detailed in Chap. 8. 


3.3 Rao-Blackwellised Particle Filter 


One of the challenges in implementing a particle filter is that the number of particles 
required to make a good approximation to the desired posterior pdf can grow expo- 
nentially with the dimension of the state space. In some circumstances, it is possible 
to mitigate this by incorporating an analytic representation of the distribution of part 
of the state given a sample of the remainder of the state. For example, suppose that 
the measurement function can be decomposed into two parts 


a = hy, (xi) + hz (x) + wi. (3.12) 
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where x! and x? are disjoint sub-vectors of the state x+. In this case we can write 
1.2 1 2,1 
p (xi. Xj|Z;, Z;-1) =p (xilz., Zi) p (xix; Zr, Z,-1) . (3.13) 


The two densities above can be estimated using different filters. When the function 
h? (x) is linear and the noise is Gaussian, the second density p (x?IxL. Zi, Zi) 
can be estimated using a Kalman filter, even if the first function hj (xi) is nonlinear. 
The state vector that needs to be sampled is then x} not [x}, x7] and the sampling 
process can use fewer samples for a given degree of accuracy. 

When a particle filter is used for the nonlinear part of the measurement problem, 
the conditioning of the second state density p (xi, Zi, Z1) leads to a separate 
Kalman filter for each particle. Each Kalman filter uses the sampled value of 
the sub-state x} as though it were the truth. This arrangement is referred to as a 
Rao-Blackwellised particle filter [15, 29, 38]. 
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Chapter 4 
Aircraft Prior Based on Primary 
Radar Data 


The Bayesian approach described in the previous chapter is a recursive method that 
calculates the posterior state distribution at each measurement time from a distri- 
bution at the previous measurement time. It fundamentally requires knowledge of 
three probability density functions: the prior distribution of the state at initialisation, 
p (x(0)); the state evolution p (x,|x,—1); and the measurement likelihood p (z;|X;.). 
Chapter 5 addresses the measurement probability density and Chaps. 6 and 7 discuss 
the state transition model. This chapter discusses the prior state distribution and the 
method used to define it. Intuitively, one might expect this prior to have a significant 
influence on the probability distribution at later times: a larger uncertainty in the 
prior might be expected to lead to a greater spread of uncertainty in the final pdf 
compared to a prior pdf with smaller uncertainty. 

In the MH370 search, there are two data sources that are available to construct 
the prior. The aircraft reports its own location and other information to the ground 
via a satellite link using the Aircraft Communications Addressing and Reporting 
System (ACARS). Data from ACARS is available for the accident flight only up to 
the point where communications were lost: the final ACARS report was at 17:07:29. 
In Chap. 9, other flights with known aircraft locations are used to validate the models 
used for the accident flight. For these, ACARS data is available and this data is used 
to construct the prior. 

The second source of prior information is radar. For the validation flights this radar 
data is not available and nor is it required given the presence of ACARS logs. For the 
accident flight, primary radar data provided by Malaysia is available from after the 
loss of communications up until 18:22:12. The radar data contains regular estimates 
of latitude, longitude and altitude at 10 s intervals from 16:42:27 to 18:01:49. A single 
additional latitude and longitude position was reported at 18:22:12. Figure 4.1 shows 
the radar data overlaid on a map. Under radar coverage, the aircraft turned sharply 
at approximately 17:24, crossed over Malaysia, and then turned to the North-West 
at 17:53. 

The Bayes filter requires a prior over the full state space, not merely latitude and 
longitude. The development to follow will lead to a state vector containing several 
other parameters. Where possible itis preferable to specify a prior on these parameters 
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Fig. 4.1 Primary and secondary radar data available for MH370 


using radar data or ACARS rather than subjective belief. Where this is not possible, 
the philosophy has been to use priors that are diffuse to avoid prejudicing the filter 
output. It is possible to derive the angle and speed of the ground velocity from the 
radar reports by assuming a simplified almost constant velocity model and applying 
a Kalman filter. This assumption is acceptable for the primary radar data because the 
reports are closely spaced in time. Figure 4.2 shows the derived speed and direction 
obtained from this filter.! The speed estimates vary dramatically during the first turn, 
which is not an accurate representation of the aircraft speed at this time. It is likely 
due to the mismatch between the assumed linear Kalman filter model and the high 
acceleration manoeuvre performed by the aircraft. Since these artefacts are localised 
to the time ofthe turn the influence on the state at the end of the sequence is negligible. 

The final reported position from radar was at very long range from the sensor and 
there was a long time delay between it and the penultimate radar report. This report is 
at long range and it is likely to have rather poor accuracy because the angular errors 
translate to large location errors at long range. The radar report at 18:22 is closer to 
the penultimate report at 18:02 than the filter speed predicts. Also, it was observed 
that the range ring derived from the timing measurements at 18:25 and 18:28 are 
closer to the 18:02 report than predictions based on either the 18:02 filtered speed or 
the 18:22 filtered speed. Figure 4.3 shows the relative positions of the 18:25 arc and 
the filter predictions based on data up to 18:02. Collectively these data points suggest 
that the aircraft may have slowed down at some point between 18:02 and 18:22. In 


l The measurement error was assumed to have a standard deviation of 0.5 nm and the process noise 


variance was 3.5 x 10^ nm? s73. The process noise was adjusted to minimise the mean squared 


prediction error. 
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Fig. 4.2. Smoothed estimates of speed and direction derived from radar data. Dotted lines show 
covariance of estimates, illustrated as mean plus and minus three-sigma value 


addition, the ground speed observed by the radar prior to 18:02 is relatively high and 
implies that the aircraft would have been at low altitude. This is likely to result in 
poor fuel efficiency, and in order to maintain flight for the duration indicated by the 
satellite data, the aircraft would have had to slow and increase altitude at some stage 
to conserve fuel. This is also consistent with a potential speed change between 18:02 
and 18:22. 

The 18:22 radar observation was not used quantitatively because the latitude and 
longitude derived from it are likely to be less accurate at long range and the aircraft 
may have manoeuvred prior to 18:22. The radar observation was deemed to indicate 
that the aircraft did not turn between 18:02 and 18:22, but the numerical values were 
not used. Instead, a prior was defined at 18:01 at the penultimate radar point using the 
output of the Kalman filter described above. The position standard deviations were 
set to 0.5 nm and the direction standard deviation to 1°. Figure 4.3 shows predictions 
ofthe mean of this prior from 18:02 to 18:25, shown in yellow, and one-sigma lines at 
+1°. The 18:22 radar point, at the end of the radar track, is clearly within the azimuth 
fan. As described above, the filtered speed at the output of the Kalman filter is not 
consistent with the 18:25 measurement and predictions based purely on this will have 
a likelihood very close to zero. In addition, the model discussed in Chap. 6 specifies 
air speed in terms of Mach number. The manoeuvre model described in Chap. 7 allows 
for speed changes and these will be randomly sampled by the proposal distribution. 
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Fig. 4.3 18:02 prediction to 18:25, shown in yellow. The Malaysian military radar track is shown 
in white, on the right. The near-vertical white line on the left corresponds to the 18:25 BTO arc 


Rather than trying to specify when the speed change occurred, the filter was expected 
to learn this information. This provides a richer description of the trajectory since the 
timing of a speed change and the new speed selected are coupled together to arrive at 
an appropriate position at 18:25. As will be seen, the filter had no difficulty finding 
paths that agreed with the measurement data. An initial Mach number was selected 
from a uniform prior between 0.73 and 0.84; this was chosen on the basis of expert 
advice to ensure that the required flight endurance is achievable. 

Itis possible to define a much earlier prior using only the ACARS data from early 
in the flight. In this case, the turns illustrated in the radar data become part of the 
unknown aircraft flight path to be estimated. Chapter 10 illustrates that this approach 
leads to a larger search zone, including the search zone resulting from the use of the 
radar data, and extending further North. This broadening occurs because the flight 
path from the final ACARS report to 18:25 is much less constrained. 
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Chapter 5 
Measurement Model, Satellite 
Communications 


The Bayesian filter discussed in Chap.3 relies on knowledge of three probability 
density functions: the state prior distribution, the state stochastic model, and the 
measurement conditional probability density. The prior used for the analysis in this 
book was discussed in Chap. 4. This chapter addresses the measurement probability 
density and Chaps. 6 and 7 discuss the state dynamics model. 

The most general measurement model was defined in (3.5) and simply states 
that the measurement is some function of the aircraft state and measurement noise. 
In many systems it is not too restrictive to assume that the noise is additive, in which 
case (3.5) becomes 

Zk = hy (xx) + Wi. (5.1) 


Prescribing the nonlinear function h; (x4) and the statistical distribution of the mea- 
surement noise provides a complete description of the measurement probability den- 
sity. The measurements available for the accident flight are timing and frequency 
logs of communication messages between the aircraft and a ground station. Details 
of the communication system software and hardware combine with the physics of the 
communication geometry to determine the nonlinear measurement function. The sta- 
tistics of the noise were determined empirically by analysing a population of real 
measurements for known aircraft states. 

This chapter gives a brief overview of the satellite communication system and then 
describes the nonlinear measurement functions and empirical noise models for the 
timing and frequency measurements. The majority of the communications messages 
available were automated signalling messages, but there were also two telephone 
calls made to the plane that remained unanswered. The first of these is particularly 
important because of when it occurred. The chapter concludes with a description of 
the measurement model for telephony. Further details may be found in [2]. 
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5.1 Satellite Communications System 


The accident aircraft was fitted with a satellite communications terminal that used 
the Inmarsat Classic Aero system [2]. This system uses a satellite to relay messages 
between the aircraft and a ground station. In the accident flight the messages were 
passed between the aircraft and a ground receiving unit located in Perth, Australia, via 
the Inmarsat-3F1 satellite. Figure 5.1 illustrates the satellite communication system 
in use during the accident flight. The aircraft is referred to as the Aircraft Earth Station 
(AES) and the ground receiving unit is referred to as the Ground Earth Station (GES). 
Inmarsat-3F 1 is a satellite in geosynchronous orbit with longitude 64.5? East and was 
used exclusively for the duration of the accident flight. 

An AES is equipped with a satellite data unit that comprises a satellite modem 
and auxiliary hardware and software. Transmission of data over the satellite is via 
bursts which are scheduled to arrive at the GES at a specified time and frequency. 
Communications from multiple users are coordinated by the allocation of different 
time and frequency slots to each user. Return channel (AES to GES) time slot bound- 
aries are referenced to the forward channel (GES to AES) [2]. The duration of each 
time slot is sufficient to account for all possible positions of the aircraft with respect 
to the satellite. The width of each frequency slot is determined by the data rate and a 
guard width that accounts for possible variations in the satellite oscillator frequency 
and other possible frequency offsets. Frequency compensations applied onboard the 
aircraft (aircraft induced Doppler pre-compensation) and at the ground station serve 
to reduce the possible difference between the expected and actual frequency of the 
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Fig. 5.1 System model of the satellite communication system 
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messages received from the aircraft. The on-ground compensation makes use of a 
second ground station located in Burum, Netherlands that transmits a reference sig- 
nal to the Inmarsat-3F1 satellite which is relayed to the Perth GES. Its purpose is 
to enable the receive modem in the GES to compensate for the Doppler frequency 
shift from the satellite to the Perth GES. This compensation process is referred to as 
Enhanced Automatic Frequency Correction. 

After the Enhanced Automatic Frequency Correction process, the expected time 
of arrival of each communications burst is compared with the actual time of arrival 
and the difference between the two is referred to as the Burst Timing Offset (BTO). 
The BTO is minimised when the elevation angle to the satellite is 90? and increases as 
the aircraft moves away from the sub-satellite position. Hence, the BTO is a measure 
of how far the aircraft is from the sub-satellite position. Similarly, the difference 
between the expected frequency of each communications burst and the actual received 
frequency is referred to as the Burst Frequency Offset (BFO). The BFO is a function of 
the relative speed between the aircraft and the satellite. Given that the satellite position 
and speed are known, the BFO provides information about the aircraft velocity vector. 
The BTO and BFO are logged by the ground station for every communications burst. 
This logging was a relatively recent addition to the ground station following the 
Air France 447 accident [2, 45] and was intended to assist in locating an aircraft. 
Statistical models for these two measurement functions are now developed. 


5.20 Burst Timing Offset 


The Inmarsat Classic Aero system allocates a time slot for communications based 
on a nominal propagation delay that assumes a nominal satellite position and a 
nominal aircraft position. The nominal aircraft position is at zero altitude directly 
below the satellite's nominal orbital position of 64.5?E longitude, zero latitude and 
an altitude of 35788.122km. The round trip delay is proportional to the distance 
from the ground station to the actual aircraft location via the actual satellite location. 
The actual satellite position is different from the nominal satellite position because 
Inmarsat-3F 1 is not exactly motionless, but rather moves in a known way in a region 
about its nominal location. The Burst Timing Offset is the additional delay after the 
start of the allocated time slot at which the message is received [2]. The BTO is thus 
the difference between the round trip message delay and the nominal delay used 
for scheduling. In addition to the propagation delay the message delay includes the 
latency of the satellite data processing unit. Denoting a BTO measurement at time f; 


as zBTO, the BTO measurement function is 


zBTO — BT (xp, sp) + RTO (5.2) 


hi S) oT Gps) = TT Ee, (5.3) 
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where 


e T (Xx, Sx) is the round trip propagation delay from the ground station to the aircraft 
via the satellite; 

e S, is the state of the satellite at time f;, that is its position and velocity in three 
dimensions, along with the satellite oscillator's state; 

e T°" is the nominal round trip delay; 

e T^?^el is a channel dependent bias term due to processing in the satellite data 
unit; 


anomaly . . F 
T, Y is an anomaly correction term discussed below; 
e wBTO is a zero mean scalar noise process with statistics to be determined from 


measurement data logs. 


The function hBTO (Xz, Sk) is assumed to be deterministic, so the measurement vari- 
ance is the variance of the noise term wBTO, The round trip delay can be expressed 
in terms of the distances from the satellite to the ground station and the aircraft as 


2 
T (Xk, Sk) = -( |H's; — g| + |H’s, — H*x;| i (5.4) 


where c is the speed of light; g is the location of the ground station; H? selects the 
position elements of the satellite state; H* selects the position elements of the aircraft 
state; and the distance |- | is the three dimensional Cartesian distance, i.e. the /? norm. 

Combining the nominal locations of the satellite and aircraft with the known 
location of the Perth ground station gives a value of T"°™ = 499,962 us. 

There are a number of different channel types used that carry different traffic types 
and have different baud rates. Communications from the aircraft to the ground are 
typically over the R- and T-data channels with C-channel used for voice telephone 
calls. Communications from the ground to the aircraft are over the P-channel. The 
channel dependent calibration term T *'?"*! is assumed to be constant over a single 
flight but can vary between flights. A fixed value for each flight assessed was empir- 
ically derived by comparing the communications logs with known aircraft positions: 
the calculated value of T°"""*! was the mean difference between the measured BTO 
and the expected BTO calculated using the known aircraft location. For the accident 
flight this calibration is only available for the time when the plane was at the tarmac 
and for the first half hour of flight. As such, values from the previous flight were also 
used in the calculation of T^^?^ne! The majority of the messages available from the 
accident flight are R1200 messages for which T°#""¢! — —4 283 us. 

The anomaly correction term ono was empirically derived through analysis 
of a collection of communications logs. For some communication messages, typi- 
cally during initial log-on, there was a very large difference between the measured 
BTO and the nominal delay. Analysis showed that rather than simple outliers, these 
anomalous BTO measurements could be corrected by a factor of N x 7,820 us 
where N is a positive integer. The origin of these anomalous BTO measurements 
has not been fully determined, but the empirical correction time is quite close to the 
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transmission interrupt clock period of 7,812.5 us and the BTO collection process 
contains quantisation. 

The channel dependent calibration and anomaly correction terms result in a resid- 
ual measurement error that is approximately Gaussian. For the R1200 messages, the 
empirically derived standard deviation of the measurement noise wel? is 29 us, and 
for R600 messages, 62 us. For anomalous R1200 messages a standard deviation of 
43 ws was used. Figure 5.2 shows a histogram of the residual BTO measurement 
errors for R1200 messages referenced to the 7 March 2014 Thamel value, and the 
moment-matched Gaussian approximation. The data used to construct the histogram 
and the empirical parameters were obtained from logs of the 20 flights of 9M-MRO 
prior to the accident flight. 

The histogram in Fig. 5.2 has an underlying mean of 10 us. This is due to the 
channel dependent calibration term Tanne! not being stationary. Over the span of a 
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day it appears constant but in the context of the 20 flights represented in Fig. 5.2 there 
is a slow variation. As discussed above, different values were fitted for each flight. 
Figure 5.3 shows a scatter plot of the BTO errors against time. The channel dependent 
calibration term Thamel was matched to the final flight before the accident flight, 
MH371 (and the beginning of the accident flight) and the BTO errors from MH371 
on 7 March 2014 are marked as red crosses. The variation in bias is sufficiently slow 
that assuming it is the same for the accident flight as the previous flight is satisfactory. 


5.3 Burst Frequency Offset 


The Burst Frequency Offset is a function of the Doppler shifts imparted on the com- 
munication signal due to the motion of the satellite and the aircraft. The relationship 
is more complicated than a direct Doppler calculation because the aircraft software 
contains Doppler compensation that offsets the Doppler shift due to the aircraft 
motion. Although the aircraft attempts to compensate for its own motion, it does this 
under the assumption that the communications satellite is in motionless geostationary 
orbit and it does not include the vertical component of the aircraft velocity (which is 
non-zero when it is ascending or descending) [2]. Since Inmarsat-3F1 is not exactly 
geostationary, the compensation is unable to completely remove Doppler effects. 
Empirical analysis of the BFO was conducted for the 20 flights of 9M-MRO prior to 
the accident flight. This analysis used the same Doppler correction software as the 
9M-MRO satellite data unit to determine the expected BFO given a known reported 
aircraft position and velocity and compared this with the observed measurements. 

Similar to the BTO, the BFO measurement function consists of a nonlinear func- 
tion of the aircraft and satellite states and several bias terms [2] 


zero = APFO cy, sk) + wee (5.5) 
ABFO (xy, sy) = AFL? (xg, sy) + AF SM" (sy) + 6 fL? (xi) + OFS (su) 
+5 f (sy) + 6 £P (x,, s), (5.6) 


where 


AF, id (Xz, S) is the uplink (aircraft to satellite) Doppler shift; 
AF, Senn (sx) is the downlink (satellite to ground station) Doppler shift; 
ô JE P (xx) is the frequency compensation applied by the aircraft; 


ô f~ (sz) is the variation in satellite translation frequency: the satellite uses a local 

oscillator to translate the carrier frequency of the message; 

e ô qe (sz) is the frequency compensation applied by the ground station receive 
chain; 

e ô fees (Xz, S) is a slowly varying bias due to errors in the aircraft and satellite 

oscillators and processing in the satellite data unit; 

were is a zero mean scalar noise process with statistics to be determined from 

measurement data logs. 
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This function was described in detail in [2], we review it briefly and elaborate where 
the analysis herein makes different modeling assumptions to those in [2]. Again, 
the function hBFO (Xz, Sk) is assumed to be deterministic, so the measurement vari- 
ance is the variance of the noise term wBFO, This is a less reliable assumption than 
for BTO because the bias term ó Abias (Xk, Sk) changes. To compensate for this, the 
measurement variance was inflated from the empirically derived wBFO variance. 
The uplink Doppler shift can be expressed as a product of the uplink frequency 
with the inner product of the relative velocity (normalised by the speed of light, c) 
and a unit vector along the relative displacement between the aircraft and the satellite 


FP (V5s, — V*x,)! (His; — H'xj) 


AFP (xx, Sk) = 
k (Xe, SI) |H’s, — H*x;| 


; (5.7) 


where F"'P is the uplink carrier frequency; V? selects the velocity elements of the 
satellite state; and V* selects the velocity elements of the aircraft state. Similarly the 
downlink Doppler shift and the frequency compensation are given by 


p down (Vs)! (H's; — g) 


AFE (sz) = Bs a ^ (5.8) 
up (Vx T H^ E 
5 "(my = PLU the E : (5.9) 
ye 


where § is the nominal satellite position assumed by the aircraft, and F°" is down- 
link carrier frequency. The aircraft frequency compensation term ó fe (xx) is deter- 
mined using the aircraft's own knowledge of its position and velocity but with an 
assumed altitude of zero and an assumed vertical speed of zero. The modified posi- 
tion matrix H* selects only the horizontal location variables and sets the altitude to 
zero, and similarly for V*. The compensation also assumes a motionless satellite at 
its nominal satellite location of 64.5?E. The satellite altitude used in the correction 
is 422 km higher than the nominal 35788.12 km value. 

The satellite translates the frequency of messages using a local oscillator that 
is maintained in a temperature-controlled enclosure to improve its stability. During 
eclipse periods, when the satellite passes through the Earth's shadow, the satellite 
temperature drops, resulting in a small variation in translation frequency [2]. An 
eclipse period occurred during the accident flight and some of the validation flights 
were also affected by eclipses. The oscillator temperature also varies with time of 
day as the satellite orientation to the sun changes and as the temperature control 
system applies its controls. All of these thermal effects are included in the term 
6 f£? (sp). The specific details of the functions that define 6 f$% (s) and 6 f^ FC (s,) 
are proprietary of Inmarsat. 

The bias term ó tes (Xz, S) is time varying. In the BTO case the variations in bias 
were slow enough to be ignored within a single flight and we were able to assume 
that 7°h4""e! was constant for each flight. This is not the case for the BFO bias term. 
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Fig. 5.4 Results for the 2-Mar-2014 flight from Mumbai to Kuala Lumpur 


The mean bias is different between flights and even within a single flight there is 
evidence of structured variation. Figure 5.4 shows an example of the BFO measure- 
ment errors for a flight between Mumbai (BOM) and Kuala Lumpur (KUL) on 2 
March 2014. The figure shows the difference between the measured BFO values and 
predicted values (using the actual SDU software for determining ó gom (xx) in the 
SATCOM system model) based on the known geometry and aircraft velocity vector 
at the time. The bias used for the plot was obtained by analysing BFO measurements 
while the aircraft was on the tarmac. The residual error is clearly not zero-mean, and 
the mean varies with time. Substantial effort was made to characterise this structured 
bias. It was found to have a geographic dependency but it has not been possible to 
determine a quantitative function to compensate for this change in bias. 

The variations in bias shown in Fig. 5.4 happen over a timescale of minutes rather 
than hours. In the accident flight the available BFO values are generally at least an 
hour apart. This is a relatively long time compared with the correlation structure of 
the error, so the model does not use a coloured noise model for the BFO. However, 
the drift of the BFO bias means that it is not sufficient to assume that 6 feas (Xk, Sk) 
will be the same in flight as on the tarmac before takeoff. The potential variations 
were incorporated by modeling the BFO bias as an unknown constant with a prior 
mean given by the tarmac value and a standard deviation of 25 Hz. Since the BFO 
measurement Eq. (5.6) is linear in the bias its distribution conditioned on the other 
states can be estimated with a Kalman filter. This is the Rao-Blackwellised particle 
filter described in Sect. 3.3. 

Empirical statistics of the residual measurement noise wBFO were determined 
using the previous 20 flights of 9M-MRO. Data points corresponding to when the 
aircraft was climbing or descending were excluded. Table5.1 shows the empirical 
statistics of the BFO measurements for R1200 and R600 messages. The ‘in-flight 
only' statistics show the combined effects of noise and bias variation without the 
influence of ‘on-tarmac’ outliers (potentially due to taxiing). The ‘including tarmac’ 
statistics on the other hand are also influenced by the BFO bias value applied to 
keep the BFO error at the source tarmac for R1200 messages close to zero. The 
mean BFO error was close to zero in all cases, indicating appropriate 6 Jen (Xx, Sx) 
values were chosen for each flight. The statistics show that even when outliers are 
discarded a standard deviation of about 4.3 Hz is applicable. As discussed above, 
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Table 5.1 Statistics of BFO errors for 20 flights of aircraft 9M-MRO prior to MH370 


Mean BFO error | Standard Mean BFO error | Standard 
(outliers deviation of BFO | (outliers deviation of BFO 
included) (Hz) (outliers excluded) (Hz) (outliers 
included) (Hz) excluded) (Hz) 
Including tarmac | 0.2246 4.9592 0.2745 4.0192 
In-flight only 0.1079 5.4840 0.1755 4.3177 
Fig.5.5 Histogram of BFO 140 
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to be conservative and allow for potential variation in the 6 yes (Xz, Sk) value on the 
accident flight, our model assumes a noise standard deviation of 7 Hz. Section5.5 
illustrates the sensitivity of the BTO and BFO measurements to variations in the 
aircraft state. 

Figure 5.5 shows a histogram of the 3392 in-flight BFO errors. On-tarmac BFO 
errors were excluded due to the pre-biasing described above. A Gaussian fit to the 
distribution is shown as a black line. It can be seen that the distribution shows 
some non-Gaussian features and the tails of the distribution for negative errors are 
somewhat heavier than those for positive errors. 


5.4 C-Channel Telephone Calls 


There were two unanswered telephone calls from the ground to MH370 after the loss 
of radar data. These communications use the C-channel and result in measurements 
of BFO but not BTO. Initially the C-channel data was not included in the flight 
prediction but analysis from DST Group highlighted that the first of these calls 
provides critical information. The first call occurred from 18:39:53 to 18:40:56 and 
is important because the measured BFO is significantly different from the BFO on 
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Fig. 5.6 Implied track angles for MH370 with different assumed ground speeds at 18:40, 7 March 
2014 


the R1200 measurement preceding it at 18:28:15. The R1200 BFO value is consistent 
with the speed and direction of the aircraft while under radar coverage whereas the 
later C-channel BFO value is not. Assuming that the change in BFO implies a turn, 
the difference between the BFO predicted by using a MATLAB model of the SDU 
software! and the measured BFO on the C-channel was analysed as a function of 
post-turn direction and for a range of aircraft speeds and turn times between 18:28:15 
and 18:39:53. Figure 5.6 shows the residual error and it clearly demonstrates that only 
Southerly track angles are consistent with the C-channel measurements. The model 
predicted BFO values of Northerly paths are more than 10 standard deviations away 
from the measured BFO. 

The BTO measurements from the R1200 messages at 18:28:15 and 19:41:03 are 
not consistent with the velocity vector before 18:28:15. The only way to satisfy 
these measurements and maintain a feasible air speed is for the aircraft to have 
turned. However, the time window for this turn is more than an hour. The 18:39:53 
C-channel measurement is critical because when combined with the 18:28:15 BFO 
measurement it significantly restricts this turning time window to a little over 11 
minutes. Using the C-channel data restricts the aircraft trajectories much more tightly 
than using only R-channel data. 


! Note The difference between the MATLAB model output and the SDU software output was found 


to be inconsequential to this analysis for determining 6 Pu (xz) in the SATCOM system model. 
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5.5 Information Content of Measurements 


The information content of the BTO and BFO measurements is illustrated in Fig. 5.7. 
The figure shows a small part of the likelihood function of the BTO and BFO mea- 
surements at 19:41:02. The plots were created assuming an altitude of 30,000 ft. 
The BFO diagram used an assumed aircraft position of 1°S and 93.6°E and a bias 
of 150 Hz, which is the tarmac value for the accident flight. The BFO contour shape 
varies slowly with aircraft position. The figure used the measurement error model 
described earlier in this chapter, namely zero-mean Gaussian noises with standard 
deviation 29 us for BTO, and 7 Hz for BFO. 


BTO Measurement Likelihood (19:41) x 103 
—0.5 T TET: T 
| 12 
| 10 
E) 
g | 8 
3 
er | 
5 6 
4 
4 
-1.5 1 WL 1 1 
92:5 93 93.5 94 94.5 
Longitude (deg) 
BFO Measurement Likelihood (19:41) 

500 T iN 7, T T 

480 - | J 0.05 

460 F 1 

440 L 4 0.04 
& 420+ | 
E 400 - | 0.03 
a. 
v 380+ 1 

0.02 

360 - J 

il | 0.01 

320 } E 

0 60 120 180 240 300 360 
Track Angle (deg) 


Fig. 5.7 Likelihood functions for BTO and BFO measurements at 19:41:02 
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The diagrams show that the BTO provides reasonable localisation along a circle 
of a given range from the satellite. The information provided by the BFO is less pre- 
cise, providing information on speed, with standard deviation on the order of 50 kn, 
and direction on the order of 40°. 
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Chapter 6 
Aircraft Cruise Dynamics 


The final piece of the Bayesian filter is the dynamics model. In object tracking, it is 
common to model either the velocity or the acceleration of the object as a random 
walk in two or three dimensions (e.g., [7]). Such a model has the advantage of 
mathematical simplicity and can be adequate, particularly in applications where the 
delay between measurements is low and the area of interest is small. However, an 
integrated random walk is not adequate for prediction over long time periods, which 
is required for applications where the delay between measurements is significant. In 
the case of MH370, the area of interest is vast and the measurements are extremely 
sparse in time, occurring approximately hourly. The aircraft is believed to have turned 
at least once after loss of radar contact, so the dynamics model needs to be adaptive 
enough to accommodate realistic aircraft manoeuvres and yet not so loose as to 
cause unrealistic spreading of the state posterior pdf over the whole BTO arc. This 
is not possible with simple integrated random walk models, so more sophisticated 
approaches are required. 

The dynamics model that was adopted for this analysis consists of a random 
sequence of deliberate manoeuvres interspersed with periods of cruise, in which 
the speed and direction are almost constant. This model is motivated by [12, 16]. 
Commercial aircraft typically operate by setting a desired air speed in terms of a 
Mach number, which is the ratio of the aircraft speed relative to the air mass around 
it to the speed of sound. During cruise the aircraft autopilot system seeks to maintain 
a constant Mach number, altitude and control angle (e.g., heading or track), though 
random fluctuations will occur in each of these. 

Three types of manoeuvre were considered: coordinated turns, where the aircraft 
maintains a constant speed and changes control angle; altitude changes that occur 
under constant control angle and speed; and accelerations, where the aircraft main- 
tains a constant control angle and altitude and changes speed (Mach number). Details 
of the statistical models used for manoeuvres are presented in Chap. 7, while details 
about the cruise model are provided in this chapter. 
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6.1 Ornstein-Uhlenbeck Process 


As discussed above, it would be inappropriate to assume a random walk in either 
speed or direction because, over the long duration between measurements, such a 
model would either apply significant probability to speeds which are kinematically 
infeasible (if the noise strength was significant), or fail to model the statistical vari- 
ability that flight paths tend to exhibit (if the noise strength was very small). Instead an 
Ornstein-Uhlenbeck (OU) process [41] was adopted to model the speed and direction 
under control. Under this model, random perturbations are permitted, but the values 
tend to revert back to a prescribed set point. Denoting the quantity being controlled 
as x(t) and the prescribed set point as x(t), the model is specified in continuous 
time as 

dx(t) = —B [x(t) — x(t)] dt + dx (t) + dv(r). (6.1) 


The parameters of the model are the reversion rate £; the strength q of the Brownian 
motion process v(t); and the set point x (t). The set point is accommodated through the 
sampling procedure described in Chap. 7, while the reversion rate and noise strength 
are selected using a maximum likelihood procedure operating on logged aircraft flight 
data from other flights using the same aircraft type, as will be described shortly. 
For a linear, time invariant, continuous time stochastic differential equation of the 
form 
dx(t) = Ax(t)dt + dv(t), (6.2) 


the equivalent discrete process is 
Xk = Ó,x, 1 + Vi, (6.3) 
where 6, is the system matrix, calculated using the matrix exponential 
9, = explA (tk — tk-1)], (6.4) 


and the covariance of the noise v; is 
tk 
Q = | exp[A (tt — T)lq(T)(expPĪA (tk — c)]' dr, (6.5) 
th-1 


where q(T) is again the strength of the vector Brownian motion [30, 41], here assumed 
constant. In this general case x is a vector and 6 is a matrix. For the OU process, we 
consider a scalar state x and a scalar system matrix ® 


Xy = Xy + Dy Q1 — Xia) + wi, (6.6) 


where x, = X(t) and 
Pr = exp[— B(ty — te-1)]. (6.7) 
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The variance of the process noise w; is: 


Or = 35 exp[—2B(th — tk-1)]). (6.8) 


If the difference in time between two samples is large, then the new value of the OU 
process is only weakly correlated with the previous value. This can be seen by taking 
the limit of (6.7) and (6.8) as the difference in time increases: the expected value 
converges to xg, and the variance converges to a steady state value of Ax. 


When an OU process is used to model velocity, the position that results from inte- 
grating the velocity is an Integrated Ornstein-Uhlenbeck process. Analytic expres- 
sions for the parameters of this process can also be determined in closed form. In 
the present context, the overall velocity is constructed as a nonlinear composition 
of several OU processes, as described below. Consequently, position is determined 
through numerical integration of the velocity over small time steps. 


6.1.1 Determining Process Parameters 


Logged aircraft flight data from several flights were analysed in order to select the 
parameters £ and q of the OU processes for Mach number, wind, and angle. The log 
data includes variables such as position, velocity, Mach number, heading, and local 
wind velocity. In particular the log contains each quantity that is modelled using the 
OU process, namely Mach number, wind, and direction. This means that a maximum 
likelihood process can be applied directly (e.g., [37, 44]), rather than requiring more 
complicated system identification techniques (e.g., [27]). 

In the first step, the data was manually segmented into trajectory parts that do 
not include manoeuvres (such as altitude, speed or direction changes). This yields 
L sequences, where the /th sequence is (x;,1, ..., X1,K,), and the time of element 
Xj, 18 tık. The unknown parameters are the reversion rate f, the noise strength q 
and the nominal set point x; for each sequence (which is treated as a non-random 
nuisance parameter). Defining 0 = (£, q, X1, ..., Xr), the overall log likelihood of 
the observed data is: 


L Kı 
s(@) = 2. L DGxka) + » og pialtu] (6.9) 
l=1 k=2 
TN: ae nnm 
idm SED fioserow 4 Era = Prarie — A= 9,08] | 
I=1 k=2 Qux 


(6.10) 


where 4; ,, and Q; are given by (6.7) and (6.8) respectively, and C is a constant 
that does not depend on the parameter vector 0. Given a value of £, it is possible to 
maximise over (q, x1, ..., Xr) in closed form, thus parameter identification can be 
performed simply using a line search. 
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In reality, the data stored in the aircraft log files are not exact values, but rather 
filtered estimates provided by the on-board navigation system. Since the complete 
information state of the estimators is not provided, the time correlation structure is 
unknown. However, our interest in the present context is in predictions over sig- 
nificant periods of time, so log data was sub-sampled such that the time between 
samples ranges between 60 and 300s. It is assumed that, over this period, the time 
correlation of the estimates has decayed to such a point that it is dominated by the 
natural variability of the process. 


6.2 Mach Number 


In the most common modes of operation, commercial aircraft maintain an approxi- 
mately constant Mach number, i.e., the ratio between the true air speed of the aircraft, 
and the speed of sound in the airmass around the aircraft. The speed of sound in an 
ideal gas is dependent only on the air temperature: 


RT 
Csouna(T) = |J i. (6.11) 


where y ~ 1.40 is the adiabatic index, R ~ 8.314J mol ! K^! is the molar gas 
constant, M is the molar mass of the gas (~ 0.02896 kg mol`! for dry air), and T is 
the absolute temperature (in Kelvins) [26]. 

Thus the true air speed of the aircraft is given by 


Vair = M X Csound(T) (6.12) 


where T is the temperature at the aircraft location and m is the Mach number. The 
consequence of the dependence of speed of sound on temperature is that changes in 
temperature (including those due to changes in altitude) cause changes in air speed. 
The air temperature data used in this study was provided by the Australian Bureau of 
Meteorology, from the Australian Community Climate and Earth-System Simulator 
(ACCESS) Global model, referred to as ACCESS-G, outlined in [6]. Further details 
of this data source are discussed in Sect. 6.4. 

Even though the aircraft may seek to hold the Mach number constant, fluctuations 
can and do occur, e.g., due to turbulence and imperfect control. An OU model is used 
to represent this variability. The parameters Bmach = 1.058 x 107? and Qmach = 
2.05 x 1077 s^! were calculated using the procedure in Sect. 6.1.1 based on historical 
logged flight data, and represent the combination of these sources of fluctuations. 
The steady state standard deviation under this model is 3.113 x 10~*. Note that Mach 
number is dimensionless. 
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6.2.1 Cost Index 


Another commonly used mode of speed control is Cost Index. In this mode, the air 
speed is selected dynamically to optimise a combination of fuel consumption and 
travel time [34]. The pilot specifies a Cost Index value, which indicates the relative 
importance of fuel and travel time. The autopilot system then dynamically selects a 
Mach number based on the aircraft weight and altitude. This mode was simulated 
by randomly sampling a Cost Index value and altitude, and then using lookup tables 
based on proprietary data provided by Boeing. 

The effect of the dynamic optimisation is that the air speed will tend to drop over 
time if the altitude is held constant. The effect is more noticeable at lower altitudes. 
During normal operation it is typical for the pilot to initiate climbs in altitude and 
the effect is less noticeable. The effect is not evident in the validation flights studied 
in Chap. 9 because in each of these flights the aircraft altitude increases with time. 
For this reason, the majority of the experiments presented in this book do not use 
the Cost Index mode. Additional experiments were performed using the Cost Index 
mode for the accident flight. The result was a slight narrowing of the pdf as compared 
to the Mach number model, as discussed in Sect. 10.7. 


6.3 Control Angle 


The motion of an aircraft relative to the Earth is influenced by the weather in its 
environment, particularly the wind. The aircraft heading (i.e., the direction in which 
the nose points) and air speed determine its motion relative to the local air mass but 
the motion relative to the Earth is a vector sum of the air velocity and wind velocity 


Vground = Vair + Vwind: (6.13) 


as shown in Fig.6.1. We refer to the bearing angle of the ground velocity as the 
aircraft's track; elsewhere the term course is also used. It is the track that determines 
motion relative to locations on the Earth's surface, not the direction in which the 
aircraft's nose points. Similarly, the ground speed is determined from the magnitude 
of the vector sum in (6.13) and Fig. 6.1. It is generally assumed in this book that the 
aircraft can drive a control loop to influence Mach number and hence air speed but 
that the resulting ground speed is somewhat at the mercy of the wind. 

The control angle is assumed to follow an OU process as described above for the 
Mach number. The model parameters for the control angle OU process (in units of 
radians) are Bangle = 9.792 x 10-3, and dangle = 4.074 x 1078 rad? s-!. 
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Between aircraft manoeuvres, the aircraft motion is notionally straight and level. 
Commercial aircraft are generally designed to operate under autopilot assistance and 
it would be highly abnormal for a pilot to maintain direct control of this type of 
aircraft for any extended period. There are five different modes in which the aircraft 
autopilot can be programmed for steady flight. The first four of these correspond to 
different definitions of the aircraft control angle 0 (t). The last is used for navigation 
to a specified location. 

The different definitions of control angle enumerate the combinations of main- 
taining a steady heading or track relative to a magnetic or geographic bearing; these 
angles are defined in Fig. 6.1. The direction ofthe Earth's magnetic field is not aligned 
with the geographic poles and so angles measured with a magnetic compass are not 
the same as angles measured with respect to the geographic poles. The difference 
between the two angles is referred to as magnetic declination [11]. Angles measured 
with respect to the poles are referred to as true angles whereas angles measured by 
compass are referred to as magnetic angles. Magnetic declination varies with posi- 
tion on the Earth’s surface, so a particular magnetic angle corresponds to different 
true angles at different locations. Figure 6.2 shows a map of magnetic declination for 


Fig. 6.1 Ground velocity is 
determined from the vector 
sum of the velocity relative bu 
to the local air mass, and the 
velocity of the local air mass, 
i.e., the wind. The heading is 
the direction in which the 
aircraft is pointing, which is 
notionally the same as the 
direction in which it is 
moving relative to the local 
air mass, 0,. The track is the 


wind velocity 
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Fig. 6.2 Magnetic 
declination (degrees) at 10 
40,000 ft on 07 March, 2014. 
Source NOAA [31] 
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the Indian Ocean region. There are significant variations in magnetic declination in 
the areas through which the accident flight is thought to have flown. If the aircraft 
had been following a fixed compass reading then it would appear to gradually turn. 
Each of the different autopilot modes is now briefly described. 


6.3.1 Constant Magnetic Heading 


In Constant Magnetic Heading mode, the aircraft aims to hold a steady air speed 
and magnetic heading, so variations in ground velocity will occur due to variations 
in wind, and magnetic declination. The bearing of the air velocity is determined by 
adding the magnetic declination $ (x) to the control angle 


bair (t) = O(t) + 9 (x0). (6.14) 
The North and East ground velocity components follow directly from the vector sum, 


cos (Oair (1)) 


sin (Oair (£)) l Vair (t) + Vwina (t). (6.15) 


Vground (£) = | 
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6.3.2 Constant True Heading 


In Constant True Heading mode, the aircraft aims to hold a steady air velocity, 
i.e., air speed and true (non-magnetic) heading. This is the same as Constant Magnetic 
Heading except that the magnetic declination is not applied to the control angle, 


Bair t) = 0 (6). (6.16) 


6.3.3 Constant Magnetic Track 


In Constant Magnetic Track mode, the aircraft aims to hold the direction of the ground 
velocity vector, relative to magnetic North, constant. The autopilot automatically 
adjusts the aircraft heading to compensate for the wind so that the vector sum of the 
air velocity and the wind velocity aligns with the desired track angle. The aircraft 
measures the wind, so the direction and magnitude of the wind velocity vector (i.e. 
Owina and vwina) are known. The air speed vair is controlled as is the ground velocity 
angle Oground- Once again the control angle is corrected by the magnetic declination 
to give the true track angle 


Aground(t) = A(t) + é(x(r)). (6.17) 


With reference to the triangle geometry in Fig.6.1, the lengths of two of the sides 
are known and the angle between the track angle and the wind angle defines the 
angle shown at the top of the triangle. The ground speed can be determined using 
the cosine rule, resulting in: 


Vground(/) = cos [ ground (t) — Owind (t)] Vwind (t) 


+ V var (D? — vwina(t)? sin? [Bgrouna(t) — Awina(t)]- — (6.18) 


6.3.4 Constant True Track 


In Constant True Track mode, the aircraft aims to hold the true direction of the ground 
velocity vector constant. In terms of the model this is the same as Constant Magnetic 
Track except that the magnetic declination is not applied to the control angle, 


Ogrouna(t) = 0 (t). (6.19) 
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6.3.5 Lateral Navigation 


The fifth mode is referred to as lateral navigation and implies that the aircraft follows 
a great elliptical path to a prescribed destination. This is the type of navigation 
typically used between way points. In this mode a fixed latitude and longitude would 
be defined and the autopilot would determine the geodesic (shortest constant altitude 
path) linking this destination with the current location. Under lateral navigation, the 
aircraft heading and track angles will gradually change. In the present context, lateral 
navigation is similar to Constant True Track, except that the control angle (i.e., the 
track) varies as the aircraft moves. 

The path travelled under lateral navigation can be determined using Vincenty's 
reckoning algorithm [46]. The ground velocity is calculated given the local track 
angle as for Constant True Track, and Vincenty's algorithm is used to calculate both 
the new location and the updated track angle required to remain on the same great 
elliptical path (given a small time step size). 

Whereas other autopilot modes follow the same control angle indefinitely until 
the pilot modifies settings, lateral navigation permits the pilot to program a sequence 
of waypoints that can be either pre-defined, named waypoints, or manually entered 
coordinates. The autopilot constructs a route which passes through each in turn, fol- 
lowing great elliptical paths between them. One way of simulating lateral navigation 
would be to draw a random sequence of waypoints and this would be quite success- 
ful for typical commercial flights. However, custom waypoints can be entered and 
on an atypical flight could even be likely. In the present context, these deliberate 
manoeuvres are modelled as changes in the control angle, as with other modes. 

Expert advice indicates that if the autopilot system is operating in lateral nav- 
igation mode and it reaches the final programmed waypoint, then it reverts to the 
previously selected heading hold mode, i.e., Constant Magnetic Heading or Constant 
True Heading. This behaviour has been modelled by simulating a process which 
switches from lateral navigation to either Constant Magnetic Heading or Constant 
True Heading at a random time. This is effectively a special manoeuvre that can 
only happen once. The probability that the switch occurs at time t is modelled as 
exponentially distributed with a mean time of 6/In(2) hours, which means that the 
probability of switching in 6h is 0.5. The exponential model is discussed in detail in 
the next chapter. 


6.4 Wind 


As described in the preceding sections, the aircraft controls air speed and either head- 
ing ortrack angle. In Constant Magnetic/True Heading modes, the wind will influence 
both the ground speed and the ground track angle. In Constant Magnetic/True Track 
and lateral navigation modes, the wind will influence the ground speed. For this 
reason, wind must be incorporated into the analysis. 
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Fig. 6.3 Wind data for 8 March 2014, 00:00, at 175 hPa (approximately 41,000 ft). Contour lines 
are drawn every 5 kn. The zonal component of the wind is the component blowing toward the East, 
while the meridional component is the component blowing toward the North 


The wind data used in this study was provided by the Australian Bureau of Mete- 
orology from the ACCESS-G model, outlined in [6]. The model provides estimates 
of wind at air pressures from 150 to 300 hPa in 25 hPa steps and from 300 to 500 hPa 
in 50hPa steps. This corresponds to an altitude range from approximately 18,000 to 
44,000 ft. The model provides predictions every three hours at latitudes and longi- 
tudes spaced by 0.375? and 0.5625? respectively. Linear interpolation is performed 
in time, altitude, latitude and longitude. An example of the model output is shown 
in Fig. 6.3. 

The predictions provided by the climate model are averages over time and space, 
so even if they are very accurate they do not exactly provide the wind experienced 
at a particular location and time. For this reason, two independent OU processes are 
used to model the error in the Eastern and Northern components of the wind data. 
Parameters of the model were determined using the process described in Sect. 6.1.1 
based on logged aircraft flight data (1.e., using the predicted wind value as the nominal 
set point). The OU process (in units of knots) parameter values were found to be 
Bwind = 0.001087, and qwing = 0.07021 kn? s^! This corresponds to a steady state 
standard deviation of 5.684 kn. 


6.5 Altitude 


The measurements available from the Inmarsat satellite are largely insensitive to the 
aircraft altitude, but modelling of altitude is necessary due to the coupling of altitude, 
air temperature and air speed (via the Mach number), and due to the variation of wind 
with altitude. This is modeled by assuming that altitude is constant during cruise, 
but permitting deliberate changes in altitude as described in Chap. 7. 
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6.6 Putting It Together 


In summary, the process for simulating aircraft cruise dynamics follows these steps: 


1. Simulate the Mach number using the OU process: 
my = ni + Pk, mach (my i — rii) + Wi mach 


where wk, mach is a zero-mean, Gaussian random variable with variance Q mach. 
and Øz mach and Qk mach are calculated with (6.7) and (6.8) using Bmach and 


mach: 
2. Simulate the control angle using the OU process: 


Ök = 2 F Px angle (Ox—1 uni 9i) + Wk,angle 


where wz angle is a zero-mean, Gaussian random variable with variance Qk angle, 
and Øk angle and Qk, angle are calculated with (6.7) and (6.8) using Bangle and 


dangle- 
3. Look up temperature in prediction table, store as Tp. Calculate true air speed as: 


| Y RT, 
V f. c be 
k air Mk M 


4. Look up wind in prediction table, store values as Vong 


true wind using the OU process: 


east 


k wind’ Simulate the 


and v 


north |  -north . north -north north 
Vi. wind = Vi, wind + Ck, wind Vi 1, wind — V&—1,wind) + Wk, wind 


Vi ind = Vewing + Dr wind Vera wind — Viet wina) + i wind 

where DO. and wi are zero-mean, Gaussian random variables with variance 
Qt. wind; and Dr wind and Q.wind are calculated with (6.7) and (6.8) using Bwing 
and qwing. 

5. Calculate ground velocity Vgroung from air speed, control angle and wind velocity 
based on mode in use, as described in Sect. 6.3. 

6. Using ground velocity, predict forward to calculate new position. 


The lateral navigation autopilot model uses a time step of 60s. Although it is 
possible to determine great elliptical paths with a very high degree of accuracy, the 
wind is spatially varying and the model samples the residual wind error, so larger 
steps are not appropriate. The other autopilot models use numerical approximations 
to propagate the state so a shorter time step is appropriate: we used a time step of 10s. 
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Chapter 7 
Aircraft Manoeuvre Dynamics 


As discussed in the previous chapter, the dynamics model used for this analysis 
consists of a sequence of deliberate manoeuvres interspersed with periods of cruise, 
in which the speed and control angle are almost constant. This chapter provides 
details of the statistical models used for manoeuvres. These models describe both 
the frequency of the aircraft manoeuvres and how the aircraft state changes as a result 
of each manoeuvre. 

Three types of manoeuvre were incorporated into the model: turns, that result 
in changes to the aircraft control angle; accelerations, that result in changes to the 
aircraft Mach number; and vertical manoeuvres, that result in changes to the aircraft 
altitude. Each type of manoeuvre was assumed to occur independently. In practice 
it is common for multiple changes to occur together. The independent model does 
not preclude this but it does not favour it either. Without data to build a correlated 
model, independence is a pragmatic assumption. 


7.1 Manoeuvre Frequency 


Each of the aircraft manoeuvres can be described by three parameters: the start time 
of the manoeuvre, the rate of manoeuvre (which will be assumed to be constant), and 
the extent or duration of the manoeuvre. For example a turn can be defined using its 
start time, the rate of change of angle, and the angle turned through. We first describe 
the model for the start time of manoeuvres, which is the same for each manoeuvre 
type. 

The time between manoeuvres was modelled as an exponential distribution. The 
exponential distribution has a single parameter, which can be interpreted as the 
average time between events. The different types of manoeuvre are modelled with 
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a single average manoeuvre period denoted by t. The use of a single time constant 
implies an intuitive model where the aircraft can be thought to be either manoeuvring 
or not, based on the level of pilot interaction. If each manoeuvre was instantaneous, 
for example the speed changed from Mach 0.83 to Mach 0.81 in zero seconds, then the 
exponential distribution gives rise to a Poisson arrival process with an average number 
of T /1 accelerations in T seconds. In practice, the model takes a finite time to change 
from the pre-manoeuvre state to the post-manoeuvre state and the exponentially 
distributed delay to the next manoeuvre is applied after the end of the previous 
manoeuvre. The impact of this is to increase the effective average time between 
manoeuvres by the average manoeuvre duration. It also imposes a maximum number 
of manoeuvres in time T whereas the Poisson distribution has infinite support. We 
assume that each manoeuvre duration is short enough to omit it from the model 
description below. 

The probability of making a turn at time ¢ under the exponential model is 
given by 


1 t 
pít; n= ze {-<}. (7.1) 
T T 
and the probability of having no turns in the interval [0, T] is 


p(t > T; t) = exp |-=] ; (7.2) 


Using this model, and neglecting the turn duration, the prior likelihood of a sequence 
of N turns at times 0 € ti <...ty < T is 


N 
ptus; t) = I] P (talti:n-1; T) | P GN > T|t:w) 
n=1 


=r exp |-=] ; (7.3) 


where fo = 0 is not a turn time but makes the notation more convenient. 

Similar expressions hold for the accelerations and vertical manoeuvres. Since 
all three manoeuvres are assumed to occur independently, the probability of a 
sequence of turns, accelerations and vertical manoeuvres is the product of the turn 
sequence probability, the acceleration sequence probability and the vertical manoeu- 
vre sequence probability. 
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7.2 Manoeuvre Extent 


The model requires a statistical description of how the aircraft manoeuvres as well 
as when it manoeuvres. As stated above, we assume that all manoeuvres happen at 
a constant rate; for example, the heading angle could change steadily at 0.5? s^. 
For real aircraft manoeuvres, the precise value of this rate can be different but in 
the context of the flight prediction required for MH370, it is possible to assume 
fixed values without losing diversity of the sampled paths. Different rates of turn 
and longitudinal accelerations were explored in earlier versions of the model but 
these changes simply resulted in a negligible increase in position uncertainty after 
the manoeuvre. 

The aircraft turns were assumed to follow a simple bank model where the rela- 
tionship between the bank angle and the angular velocity of a turn follows 


tan (Obank) S8 
(p = e; 
y 


(7.4) 


where Qbank is the aircraft bank angle and g is the acceleration due to gravity. Tight 
turns with a high angular velocity require low speed or high bank angle. Varying 
the bank angle was found to produce only minor variations in the overall aircraft 
trajectory since turns were only a small portion of the trajectories, so a fixed bank 
angle of 15? was chosen. This is nominally a steep bank for a commercial aircraft, 
although well within its performance limits. At an air speed of 500 kn, for example, 
the turn rate would be approximately 0.6? s^! and a 90? turn would take 2.5 min to 
perform. 

The final parameter required to characterise a turn is the duration of the turn, or 
equivalently the total change in angle. The aircraft was assumed to not make turns of 
more than 180° and the change in angle was uniformly sampled from — 180° to 180°. 
Itis possible for the model to sample two or more turns in very quick succession and 
as such looping turns of more than 180? are possible, though unlikely. The likelihood 
of a quick succession of turns like this depends on the mean time r. 

The accelerations were also assumed to occur with a fixed rate of change of 
Mach. The assumed rate corresponds to a change of Mach of 0.1 in one minute. 
This nominal number was not matched to aircraft performance standards but instead 
simply acts to define a finite duration for each speed change. The Mach number after 
an acceleration was uniformly sampled from 0.73 to 0.84; this range was chosen on 
the basis of expert advice to ensure that the required flight endurance is achievable. 

The vertical manoeuvres were assumed to occur at a nominal rate of 4,000 ft per 
minute. Again this rate is not meant to model the actual typical behaviour of com- 
mercial aircraft. Rather it provides a reasonable finite time extent to the manoeuvre. 
The altitude after manoeuvre was uniformly selected in steps of 1,000ft between 
25,000 and 43,000 ft. 
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Using finite time to execute manoeuvres rather than making instantaneous changes 
allows for the possibility that a measurement could have been collected during a 
manoeuvre. The model retains the ground velocity vector but not the altitude rate, 
so in principle it could describe measurements collected part way through a turn or 
acceleration but not during climbs or descents. The vertical rate will affect a BFO 
measurement but since the model only uses a nominal vertical rate it is unlikely to 
match any actual vertical manoeuvre in detail. In practice, if a measurement was 
collected during a turn or acceleration it would be very difficult for the filter to infer 
the trajectory since it would need to model the instantaneous ground speed at the 
measurement time to match the measured BFO. For the C-channel data at 18:39 and 
23:15 there are clusters of BFO measurements and these appear to be statistically 
stationary, i.e., they do not support the premise that the aircraft is turning during 
collection. 

The three types of manoeuvre are independent and the distribution of manoeu- 
vre extents is uniform for each type, the overall prior probability of a sequence of 
turns [(fe,1, 01), . .. . (fo, No, Ov, )]. accelerations [(fa,1, v1), . .. . (fa, Na» VN,)] and alti- 
tude changes [(t),1, h1), . ... (tg, v, hw, )] is given by 


p (tton Oy v, bs Ua i s VENA J Univ Pad: ) 


— exp [-z| CNN (rsa) ^ (ran) ^. (7.5) 


where 66 = 27 rad is the span of the uniform turn angle distribution, da = 0.11 
(Mach) is the span of the uniform speed distribution, and 5h = 19 is the number of 
possible discrete altitude values. 


7.2.1 Parameter Selection 


The exponential delay model is parameterised by the average time between events, 
t. It is not clear what this value should be or even if it should be constant. In 
particular one could argue that the behaviour of the aircraft during the accident flight 
did not match typical commercial aircraft. The average delay is therefore treated 
as a constant hyperparameter that is potentially different for every hypothesised 
flight trajectory, but fixed over time.! A Jeffreys prior [24] was applied to the time 
constant; this prior is a non-informative distribution in the sense that it is invariant 
to the parameterisation used. For example, it yields the same result for the two 
parameterisations of the exponential distribution, i.e., the rate parameterisation and 


! Additional experiments were conducted in which the parameter was permitted to change at 19:41 
(well after the initial manoeuvre), but little change to the final distribution was evident. 
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the time-constant parameterisation. The Jeffreys prior is proportional to the square 
root of the determinant of the Fisher information; in the case of an exponential 
distribution, it is given by 


K(t, zd <t< 
ea ee hee. (7.6) 
0 otherwise 


where K(t, T2) = (10g) — log(ri)) 7 A support range of 0.1 « c < 10h 

between events was chosen for each manoeuvre type, spanning the range of cases 

when manoeuvres occur every few minutes, to when manoeuvres are rare in the entire 

flight. A single, common parameter was drawn for all three types of manoeuvre. 
The joint pdf of the full manoeuvre description is given by 


p(t m O1:ny}s Ua iN o Viu p Un aw s Mien, ds r) 


= p(t) exp [El (se) ^ (s) ^ (ran) ^. pm 


7.2.0 Manoeuvre Model Summary 


Figure7.1 shows a graphical representation of the aircraft dynamics model. 
Ignoring measurements (which are not illustrated), the sequence of manoeuvre times 
and extents are conditionally independent from other kinds of manoeuvres given the 


Fig. 7.1 Graphical model 
showing the conditional 
dependencies in the 
probabilistic model, where 
X = (xo, ..., Xx) represents 
the entire trajectory 
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common manoeuvre time constant. The aircraft state forms a Markov chain condi- 
tioned on the manoeuvre parameters. The measurements (not illustrated) are condi- 
tionally independent from each other given the state sequence. 


7.3 Example Realisations 


In order to illustrate the behaviour of the complete dynamics model, Fig. 7.2 shows 
example realisations of trajectories randomly sampled from the model. We emphasise 
that these do not incorporate information from the SATCOM measurements, but 
rather represent samples of the prior distribution over trajectories (e.g., p(x) in (3.3)). 
Each trajectory starts with independent samples of the manoeuvre time constant and 
initial location and velocity sampled from the accident flight prior, and then draws 
an independent set of sampled turns, and speed and altitude changes. The prior starts 
the particles at 18:02 and the figure shows each trajectory through to 00:19, the time 
of the final SATCOM message. As the figure shows, some of the trajectories exhibit 
many manoeuvres, while others do not turn or change speed at all. 

The histogram of the number of turns and speed changes hypothesised by the prior 
model is shown in Fig. 7.3. The model is in effect a mixture of Poisson distributions, 
consisting of a continuum of components ranging from very low times between 
manoeuvre (six minutes) to very high times between manoeuvre (ten hours). The 
distribution is the same for each type of manoeuvre since the prior distributions of 
time constants are the same. 
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Fig. 7.2 Example trajectories sampled from prior distribution of flight dynamics: a 10 random 
realisations; b 100 random realisations 
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Fig. 7.3 Histograms of the number of manoeuvres selected by 200,000 independent realisations 
of the model (i.e., not using any BTO/BFO measurements) 


Although it is unlikely that any given trajectory will travel for more than six hours 
with no manoeuvre, those that do all end up in roughly the same place. Figure 7.4 
shows a contour representation of the pdf of latitude and longitude at 00:19 that was 
generated using 200,000 particle draws from the prior model with no measurements 
and an isotropic kernel. The red diamond in the figure shows the mean of the prior. 
There is a clear peak corresponding to paths that make no manoeuvre. This is smeared 
around an arc corresponding to paths that turn once but maintain a steady speed and 
also radially from the start point corresponding to paths that do not turn but do change 
speed. The model samples a hyperparameter for the time between manoeuvres from 
a prior that covers 0.1 to 10h. When a particle samples a very high value for the time 
constant it is unlikely to choose to manoeuvre. Where manoeuvres were selected it 
was again unlikely to choose to change speed and direction, leading to the contours 
shown. The altitude behaviour of the model is not visualised in the figure. There 
is an approximately circular region around the initialisation point where the pdf is 
elevated. This corresponds to particles that have sampled a very small time between 
manoeuvres: these turn so much that they circle back on themselves repeatedly and 
are effectively trapped close to their starting location. Beyond this circle, the pdf is 
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Fig. 7.4 Distribution of position at 00:19 using only model predictions (i.e., not using any of the 
BTO/BFO measurements). The red diamond marks the mean on the prior 


lower behind the start point than in front of it. It is also clear that the search region 
is in a very low probability part of the pdf: this indicates that the prior is not biased 
towards selecting a particular part of the 00:19arc. The search region is inside the 
lowest probability contour used in the map. 
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Chapter 8 
Particle Filter Implementation 


Solution of the Bayesian estimation method described in Chap.3 requires one to 
recursively integrate the aircraft dynamics pdf (3.6) and multiply it by the likelihood 
(3.7). Since the measurement model is highly nonlinear and the dynamics model is 
hybrid discrete-continuous, there is no way to produce a closed form posterior distri- 
bution. An alternative is to approximate the distribution numerically. As introduced 
in Sect. 3.2, the Sample-Importance-Resample (SIR) particle filter draws random 
samples from the dynamics model and weights them according to the measurement 
likelihood. This amounts to approximating the posterior distribution as 


š 
Prlz © S wié (xe — Xf) , (8.1) 


psl 


where the w? are referred to as weights (and sum to unity) and the x? are referred 
to as particles. The convergence properties of this approximation in the limit as the 
number of particles P increases have been well studied, e.g. [14, 21]. In the SIR 
version of the particle filter, the particles are randomly generated from the dynamics 
model and the weights are 


wh «x [II p(x). (8.2) 


k'ity Sty 


A problem with sampling from the dynamics is that this can be a very diffuse 
distribution. In the MH370 case, the model allows for turns and speed and altitude 
changes, and potentially several of each can be sampled between measurements. 
The proportion of particles that sample a trajectory close to the measurements will 
be small and a very large number of samples will be required to capture the high 
probability regions. This is a well known issue for filtering in high dimensional state 
spaces. 
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Resampling is one strategy that is used to improve the number of particles follow- 
ing trajectories with relatively high likelihood. It does this in a sequential manner, 
in which at each time step unlikely particles are replaced with copies of high likely 
particles though a random sampling process. Initial approaches used these conven- 
tional techniques, but it was found to be preferable to be able to process very high 
particle counts and to adaptively increase the number of particles used until it was 
possible to identify an adequate number of likely paths, rather than processing a 
pre-specified number of particles for each time step sequentially. To achieve this par- 
ticles were propagated and weighted individually; this also reduced the size of the 
data structures required and allowed preliminary results to be extracted as the filter 
was executing. The approach adopted was a form of branching mechanism which 
repeatedly constructs full trajectories. 

The method resampled each particle separately, branching a new set of particles 
from each parent instead of resampling a fixed number of particles across all of the 
particles at a given time. The branching naturally leads to an exponential growth 
in the number of particles with time and this was mitigated by pruning extremely 
unlikely paths when their likelihood became too low. This approach is not necessarily 
computationally efficient, but in this particular application it was more important to 
broadly explore the enormous state space than to minimise computational effort. 

To motivate the approach, suppose that the distribution prior to resampling is 
approximated by (8.1), such that integrals can be approximated as 


" 
J f Gu) pGu)dx, & S wt f (xp). (8.3) 


p=1 


where it is assumed that pee, wt = 1. For each particle p, draw nz > 0 copies of 


that particle, where n? is a random variable, and for each new particle X; = x}. To 
each we apply the weight 


(8.4) 


P 


Assuming that p indexes the full set of P= » jd n? new particles, it can easily be 


shown that: 


P P P 
Dd why (#) | = DEMERS OL) = Dower RL) 5) 
p=1 


Thus resampling can also be implemented through a randomised branching proce- 
dure, recursively adapting the number of particles. This permits a form of depth-first 
search, which adaptively performs more branching when likely paths result, and 
tends to prune paths which have low probability. 
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For our experiments, we chose a procedure which branches quite aggressively 
when likely paths are discovered, and prunes extremely unlikely paths. Likely paths 
are duplicated to form n branches. This is implemented by setting 


8(nf — ñ), we > 
pan)-1 5, amm (8.6) 
w,ó(n, — 1)+ (1 — w; )ô(n,), otherwise 


for n « 1. Thus, for particles P sampled from parent particle p with w? > n, 
WP = w? /ñ, and wide branching will occur, while for particles p sampled from 


parent particle p with w? < n, 5? = 1, but most commonly the sub-tree will be 
pruned. 


8.1 BFO Bias 


The BFO measurement has a bias term that was not able to be adequately calibrated, 
as discussed in Sect. 5.3. The model treats this bias as an random variable with a 
given prior density. It is possible to sample the bias along with the aircraft states but 
a more efficient implementation is to use a Rao-Blackwellised particle filter [15, 29, 
38]. Conditioned on the other states we can write a simplified BFO measurement 
model 

zBFO — gBFO + 5 fries 4 wBFO. (8.7) 
where 2F° is constant because all of the other states, such as aircraft location and 
velocity, are known because of the conditioning. This conditional measurement equa- 
tion is clearly linear in the bias and the noise is modelled as Gaussian, so the posterior 
distribution of the bias can be determined using a Kalman filter update. 


8.2 Algorithm 


In practical terms, the algorithm proceeds by repeating the following process for 
each particle: 


1. Randomly sample an average time between manoeuvres T. 
2. Randomly sample a starting state xo (position, Mach, control angle and altitude) 
from the prior at fj =18:01:49, which is described in Chap. 4. 
. Initialise the BFO bias Kalman filter. 
4. Perform the following recursion, starting with the sample at xy and measurement 
time index k = I: 


W 


a. Draw a sample of the trajectory from x,_; to xy using the hyperparameter t 
for selection of turns, speed changes and altitude changes. 
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b. Calculate the measurement likelihood p (z;|x;,) and use it to update the tra- 
jectory weight w? = p (2x|Xx) wz. 4. 

c. Use the sampled trajectory to update the BFO bias Kalman filter. 

d. If we have reached the final measurement k = K, store the trajectory and 
weight. 

e. Otherwise, if the accumulated weight is too low, i.e., we « n then branch a 
single time with probability w? and weight w? = 1, terminating the recur- 
sion branch with probability (1 — w); otherwise, branch ñ times to process 


remaining time steps with weight W? = w? /n. 


The particle weights constructed by the method are not normalised. A normalisa- 
tion step is performed when the final set of weights at the last time point is used to 
construct the required pdf. The process in step 4a, namely sampling a trajectory, is 
critical and is realised through a finite time difference implementation given by the 
following steps: 


1. Randomly sample times to make the next turn, speed change and altitude change 
2. While the current sample time is before the next measurement time, f; 


Table 8.1 State vector elements 


State vector 


Latitude Degrees 


Longitude 


Degrees 


Mach number set point 


The OU parameter, denoted m, in Sect. 6.6 


Instanteous Mach number 


Mach, denoted m; in Sect. 6.6 


Control angle set point 


The OU parameter, denoted Ox in Sect. 6.6 


Instanteous control angle 
Altitude 


Control angle, denoted 6; in Sect. 6.6 
In discrete 1,000 feet steps 


Instantaneous wind error, north 


The difference between the wind and the tabulated value, 


denoted Denn. in Sect. 6.6 


Instantaneous wind error, east 


The difference between the wind and the tabulated value, 


east ; 
denoted Ue wind 1n Sect. 6.6 


BFO bias mean 


Mean of the Kalman filter used to estimate the BFO bias 
5 bias 


BFO bias variance 


Variance of the Kalman filter used to estimate the BFO bias 


of bias 


Mean time to manoeuvre 


T 


Autopilot control angle mode 


Used to choose between constant true/magnetic 
heading/track and lateral navigation 


Cost index value, when used 


An integer between 0 and 100 
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a. If the current sample time is the time of a manoeuvre (turn, speed change or 
altitude change), then execute the manoeuvre and sample a new time to make 


the next manoeuvre 


b. Otherwise predict ahead 10s (or to the next manoeuvre or measurement, 
whichever occurs first) 


A manoeuvre is executed by making a sequence of 1s steps. For each step the 
angle, speed or altitude is incremented and the aircraft position is predicted ahead. 
The increments continue until the new desired angle, speed or altitude is achieved. 
The procedure for state prediction under cruise dynamics is summarised in Sect. 6.6. 
The state vector used for the model is given in Table8.1. There are a large number 
of parameters involved with this model and the full description of these is provided 


in Table 8.2. 


Table 8.2. Summary of filter parameters 


Initialisation 
Latitude Gaussian s.d. 0.4 min 
Longitude Gaussian s.d. 0.4 min 
Control Mach Gaussian s.d. 0.03 
Control angle Gaussian s.d. 1? 
BFO bias Gaussian s.d. 25 Hz 
Mach deviation Gaussian s.d. 0.00311 
Angle deviation Gaussian s.d. 0.0826? 
North wind deviation Gaussian s.d. 5.68 kn 
East wind deviation Gaussian s.d. 5.68 kn 
Cruise 
Mach Reversion rate Bmach 1.06 x10? 

Noise strength gmach 2.05x 1077 s7! 
Control angle Reversion rate Bangle 9.8x 10? 

Noise strength qangle 4.07 x 107? rad? s7! 
North/East wind Reversion rate wind 1.09x 10? 

Noise strength qwing 7.02x 107? kn?s! 
Manoeuvres 
All Mean manoeuvre time 7 ~Jeffreys(0.1, 10) 
Mach Uniform new Mach 0.73-0.84 
Control angle Uniform turn angle +180° 
Altitude Uniform new altitude 25,000 to 43,000 ft 
Implementation 
Branching rate Constant n within a flight 3-10 
Likelihood threshold Constant 7 within a flight e? ore ?? 
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8.3 Assumptions 


The key assumptions used by the filter are: 


1. The radar data provides an accurate estimate of the aircraft trajectory up to 
18:01:49. If, for example, the radar track used to build the prior were actually 
from a different aircraft, the predicted pdf would be invalid. Discarding the radar 
data leads to a significant broadening of the search zone, and accident investiga- 
tors believe the radar data to be correctly associated with MH370. Chapter 10.6 
considers an alternative analysis which ignores the radar data. 

2. The measurement error characteristics are known. The pdf of BTO and BFO 
measurements, in particular the standard deviation of each, is provided to the 
algorithm as a known input. Extensive study of the statistics of these measure- 
ments has been undertaken and the distributions assumed are well characterised, 
subject to the caveats discussed in Chap. 5. Incremental changes, such as minor 
inflation of the assumed BTO variance would lead to incremental changes in the 
filter output. 

3. The aircraft cruises in one of five prescribed modes and does not change between 
them (other than a single possible change from lateral navigation to constant 
magnetic/true heading). It is possible that the whole flight was continually under 
manual control but it is highly unlikely. The use of typical autopilot modes is 
reasonable. 

4. Infinite fuel: the fuel constraints on the aircraft can be applied to the pdf after- 
wards. In the simplest case, maximum reachable ranges could be used to censor 
impossible trajectories. However, analysis of candidate trajectories has indicated 
that the majority are feasible. Broad information about the fuel consumption rate 
of the aircraft has been used to inform the range of allowable Mach numbers. 

5. The fluctuations in speed, angle and the error in wind velocity are well-modelled 
by the OU process. The parameters of the OU model were selected to model these 
quantities based on recorded data from real flights. 

6. The random turn and speed change model is rich enough to describe the real 
aircraft dynamics and the implicit preferred path for the model does not bias 
prediction. Validation results in the next chapter show that the model successfully 
produces pdfs containing the true aircraft location for the available instrumented 
flights that include air speed changes, altitude changes and angle changes. 

7. The aircraft air speed is limited to the range Mach 0.73—0.84. Fuel consumption 
becomes very inefficient at speeds higher than this and at lower speeds the aircraft 
is not able to match the measurements. In practice it is likely that the viable range 
of speeds is actually much narrower than this. 

8. In Chap. 10, the pdf of the location of the aircraft at 00:19 is combined with 
a distribution of aircraft translation during descent, to give a final search zone. 


8.3 Assumptions 61 


This distribution was developed by ATSB [5] and largely determines the width of 
the search area along the 00:19 arc. It is assumed that this distribution adequately 
models the true descent scenario. 
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Chapter 9 
Validation Experiments 


The variable rate model developed for MH370 was validated by analysing data from 
a collection of flights where the true aircraft location was known; we refer to these 
as validation flights. A total of six validation flights were used for testing. Data was 
available from a larger number of flights but the majority of these were in relatively 
short segments of less than three hours. There were only a few that maintained 
communications with the satellite Inmarsat-3F1 for longer periods and it was not 
thought productive to examine the prediction performance over time segments shorter 
than three hours. Of the six flights, four are previous flights of the accident aircraft, 
9M-MRO, and the other two are flights of different aircraft that occurred at the same 
time as the accident flight. Three of the flights are relatively short and are between 
locations inside Asia, and the other three are flights from Asia to Europe. 

The data available for the accident flight consists of mostly R1200 communication 
messages at approximately one hour intervals. In order to emulate the measurement 
information content, measurement data sets were formed by randomly sub-sampling 
R1200 communication messages from the validation flights. Ten different subsets 
were formed for each validation flight, resulting in a total of sixty validation mea- 
surement sets. Multiple sets were drawn from each flight to increase the statistical 
significance of the testing data set. They also serve to illustrate the sensitivity of 
the method to the precise measurement times and values. The measurement subsets 
were selected using a randomised process that aimed to achieve an average time 
between measurements of one hour. For the analysis we treat the measurement sub- 
sets as independent Monte Carlo trials. However there are several variables that are 
in common within the group of ten subsets of a single flight: the aircraft geometry 
is obviously the same for each subset since they are drawn from the same flight; the 
residual wind errors are the same; and the BFO is known to have a slowly varying 
bias, so there can be correlation in the BFO measurements from different subsets if 
those subsets choose measurements at similar times. Finally, some subsets may in 
fact randomly choose the same measurement as another subset. 
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In each validation flight, the true aircraft location was obtained from the Aircraft 
Communications Addressing and Reporting System (ACARS) data logs. Sections 
of the flight immediately after take-off and prior to landing were not included in 
the analysis since the aircraft dynamics are very different at these times and it is 
unlikely that sparse satellite messages would be sufficient to follow it. For the longer 
flights into Europe, the aircraft changed satellites part way through the flight so 
it was not possible to use the whole flight: these were truncated near the end of 
messaging via the Indian Ocean Region satellite. The filter was initialised using the 
true aircraft location, speed and control angle with a Gaussian random error. The 
standard deviation of the initialisation error was chosen to be the same as the prior 
for the accident flight, that is 0.4? in latitude and longitude, 1? in angle and Mach 
0.03in air speed. For every subset the posterior pdf at the final measurement was 
predicted ahead to a common time, corresponding to an exact ACARS reporting 
time. This predicted pdf is compared with the ACARS report. 

This chapter first explains the particular characteristics of each flight and presents 
an example output pdf for one of the measurement sets. This output is subjectively 
compared with the ACARS truth. The statistical analysis is then presented using 
an objective performance measure over the sixty validation subsets. Table 9.1 lists 
the six validation flights used for the analysis and gives comments on some of the 
characteristics of each. The flights are ordered by time. 


Table 9.1 Summary of validation flights 


Flight path Date 


Kuala Lumpur to | 26 February 
Amsterdam 


Duration (h:mm) 


9M-MRO 


Comments 


Eclipse 


Mumbai to Kuala 
Lumpur 


Short and almost 
straight with a 
gradual late veer, 
outlier BFO 
measurements 


Kuala Lumpur to | 6 March 4:25 Yes Single climb, 
Beijing several S-turns 
Beijing to Kuala | 7 March 4:55 Yes Significant 
Lumpur climbs, Mach 
changes and 
turns, contains 
anomalous BTO 
measurements 
Kuala Lumpur to | 7 March 7:50 No Large S-turns 
Amsterdam 
Kuala Lumpur to | 7 March 7:03 No Mid flight 
Frankfurt heading 
deviations, outlier 
BFO 


measurements 
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9.1 9M-MRO 26 February 2014 Kuala Lumpur 
to Amsterdam 


The first validation flight was from Kuala Lumpur to Amsterdam on 26 February 
2014. This flight was around 7.5h long but is relatively straight. Figure9.1 sum- 
marises the features of the flight: the upper plot shows a geographic plan; the lower 
three plots show the aircraft altitude as a function of time, the aircraft heading as a 
function of time, and the aircraft Mach number as a function of time. Vertical dotted 
lines show the start and end of the time segment selected for the test. This flight con- 
tained an eclipse event so the validation also supports the Inmarsat eclipse correction 
[2]. 

Figure 9.2 shows the filtered pdf for the Kuala Lumpur to Amsterdam flight visu- 
alised using a three dimensional representation in Google Earth. The filter pdf is 
defined over a high dimensional space but for visualisation we examine the marginal 
position distribution in latitude and longitude. Because the BTO measurement error 
is relatively small the position distribution is centred on an arc of zero BTO error 
and has a narrow off-arc width. For the visualisation we marginalise the distribu- 
tion onto the zero B TO error arc and encode the probability density for each point 
along the arc using altitude: points on the curve higher above the earth correspond 
to higher probability. A white curve on the map marks the ACARS reported aircraft 
location, a yellow marker denotes the location of the aircraft at each measurement 
time. The figure also shows a representative selection of the paths sampled by the 
filter. The selection shows the highest probability path arriving at each point around 
the arc: the colour of the path shows the marginal probability at that location on the 
arc (using a colour map similar to Fig. 5.7, i.e., blue is least likely, red is most likely). 

There are a number of paths that end in significantly different locations to the 
truth. These occur because in this flight the aircraft travels in a direction that is 
almost horizontally radial from the satellite. While the aircraft moves towards the 
satellite its initial dynamics constrain the plausible paths but once it passes through 
the point of closest approach and begins to move away then it is possible to make 
turns that result in different near-radial paths. The support of these ambiguous paths is 
disjoint because of the finite number of samples: the true underlying pdf has support 
all the way around the arc. Without dynamic constraints the location of the peak of 
the pdf is simply a function of measurement noise. 


9.2 9M-MRO 2 March 2014 Mumbai to Kuala Lumpur 


The flight from Mumbai to Kuala Lumpur is the shortest validation flight selected. 
Figure 9.3 summarises the features of the flight: there is a single minor altitude change 
and the Mach number remains relatively constant. The aircraft heading gradually 
reduces for most of the flight, turning the aircraft more to the North but a veer near 
the end turns it back to the South-East. The BFO measurements for this flight contain 
several outliers that are more than 30 Hz away from other measurements at similar 
times. 
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Fig. 9.1 Validation flight 26 February 2014 Kuala Lumpur to Amsterdam. Vertical dotted lines 
show the start and end times of the segment used for validation 
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Fig.9.2 Validation flight 26 February 2014 Kuala Lumpur to Amsterdam 
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Fig. 9.3 Validation flight 2 March 2014 Mumbai to Kuala Lumpur. Vertical dotted lines show the 
start and end times of the segment used for validation 
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Figure 9.4 shows the pdf output from the filter; the true ACARS aircraft location 
is again under the main peak of the pdf. The pdf might appear to be relatively spread 
compared with some of the other flights, but the scale is much smaller in this case 
because the flight is short. 


9.3 9M-MRO 6 March 2014 Kuala Lumpur to Beijing 


This flight is the MH370 route from Kuala Lumpur to Beijing that was flown by 
the accident aircraft 9M-MRO on 6 March 2014, i.e., the day prior to the accident 
flight. Figure 9.5 summarises the features of the flight: the flight contained a single 
altitude change and several turns. Observe that there are several times where the 
heading changes for a short time before reverting back to the previous long-term 
value. These course corrections have the impact of translating the flight path and 
then returning to the previous ground velocity vector: in effect they are a kind of 
S-turn manoeuvre. If one or more of these corrections occurs between measurements 
then the most likely paths can be biased because there are no measurements to hint 
that the manoeuvres have occurred and the S-turn trajectory is less probable under 
the dynamics model than a constant angle path. 

Figure 9.6 shows the pdf from the filter: the pdf is multi-modal with three main 
peaks that are somewhat blurred together. There was a heading change just before the 
last measurement and the lack of future data makes it impossible to resolve exactly 
what manoeuvre led to the change in range rate. One of the peaks of the pdf is clearly 
centred close to the true location. 


9.4 9M-MRO 7 March 2014 Beijing to Kuala Lumpur 


This flight is the MH371 route that was flown by the accident aircraft 9M-MRO on the 
morning of 7 March 2014 and is the return flight from Beijing back to Kuala Lumpur. 
Figure9.7 summarises the features of the flight: there were three altitude changes 
and two main heading changes, the first of which was almost immediately after the 
start of the validation segment. This flight does not contain the S-turn manoeuvres 
that were present in the previous flight. In addition to the altitude changes the Mach 
number of the aircraft changed from 0.83 to 0.82. Each of these leads to a change 
in air speed. This flight contained several anomalous BTO measurements that were 
corrected using the empirical adjustment described in Chap. 5. 

Figure 9.8 shows the pdf output from the filter; the true ACARS aircraft location 
is again under the main peak of the pdf. The peak is more spread because the altitude 
changes and Mach change modify the radial speed between the aircraft and the 
satellite. The resulting BFO measurements can also be explained by course changes: 
the aircraft could change speed or it could turn slightly. The BFO measurement is not 
informative enough to discriminate strongly between these and there is not enough 
subsequent data to see which is more consistent with BTO progression. 
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Fig.9.4 Validation flight 2 March 2014 Mumbai to Kuala Lumpur 
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Fig.9.5 Validation flight 6 March 2014 Kuala Lumpur to Beijing. Vertical dotted lines show the 
start and end times of the segment used for validation 
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Fig.9.6 Validation flight 6 March 2014 Kuala Lumpur to Beijing 
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Fig.9.8 Validation flight 7 March 2014 Beijing to Kuala Lumpur 
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9.5 7 March 2014 Kuala Lumpur to Amsterdam 


This flight was from Kuala Lumpur to Amsterdam and is the same flight path as the 
first validation flight but with a different aircraft. Figure 9.9 summarises the features 
of the flight: the aircraft climbs with a sequence of vertical manoeuvres and there is 
a large S-turn manoeuvre near to the end of the analysed flight segment. 

Figure 9.10 shows the pdf output from the filter. The true ACARS aircraft location 
is under the main peak of the pdf but in this case the true location is lower in the 
tails than in the other cases. The numerical results that follow in Sect. 9.7 show that 
this flight had the worst overall performance of the validation flights, although, as 
discussed in the next section, for each subset of measurements, the final location 
is within the region containing 85 % of the probability distribution, i.e., the highest 
posterior density (HPD) interval, discussed further in Sect. 9.7.2. 


9.6 7 March 2014 Kuala Lumpur to Frankfurt 


The final validation flight was from Kuala Lumpur to Frankfurt. Figure9.11 sum- 
marises the features of the flight. It shows the full flight path, but the communica- 
tions satellite changes part way through and the test section finishes where the box 
is marked on the map. No Mach information was available for this flight. There 
was a large heading deviation mid-flight, but the aircraft eventually reverted back 
to the earlier heading: this kind of compound manoeuvre is difficult for the filter to 
characterise. This flight also contained outlier BFO measurements. 

Figure 9.12 shows the pdf output from the filter. The performance on this flight is 
quite similar to the Kuala Lumpur to Amsterdam flights. The filter has again identified 
ambiguous paths due to the relative geometry. 


9.7 Quantitative Analysis 


The examples above present a qualitative measure of performance but a more rig- 
orous objective measure is required to provide a statistical assessment of the filter 
output. So far we have been satisfied that the true location has been in an area of 
reasonable support for the pdf, but is the spread of the pdf appropriate and is the 
mode of the distribution biased? Answers to questions such as these require a much 
larger ensemble of test data. However, it has not been feasible to collect the required 
test measurements for dozens of different flights. In order to increase our confidence 
in the performance for the relatively small set of flights that is available, multi- 
ple communication measurement sets were generated for each flight by randomly 
selecting individual R1200 messages from the communication logs of each flight. 


76 9 Validation Experiments 


latitude (degrees) 


0 l 1 l l l 1 l 
30 40 50 60 70 80 90 100 110 120 


longitude (degrees) 


BR 
© 
T 
l 


w 
© 
T 
1 


altitude (kfeet) 


N 
© 
T 
I 


1 
18:00 21:00 00:00 


heading (degrees true) 


1 
18:00 21:00 00:00 


o 
e 


Mach number 
o 
oo 
T 


1 L L 
18:00 21:00 00:00 
time 


o 
a 


Fig. 9.9 Validation flight 7 March 2014 Kuala Lumpur to Amsterdam. Vertical dotted lines show 
the start and end times of the segment used for validation 
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Fig.9.12 Validation flight 7 March 2014 Kuala Lumpur to Frankfurt 
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The selection process was repeated 10 times for each flight and these 10 measure- 
ment sets are treated as independent Monte Carlo random trials for a fixed true air- 
craft trajectory. As discussed in Chap. 5, the BFO measurement errors are not truly 
independent over short time periods, which somewhat compromises the assumed 
independence. However, the common geometry of multiple sets from a single flight 
is the dominant source of correlation amongst single-flight predictions. 

We now briefly review the method used to select individual messages and the 
performance measure used for this analysis. The chapter concludes with numerical 
results from these sixty measurement sets. 


9.7.1 Measurement Selection 


The start and end time for analysis was manually selected for each flight. These times 
were chosen to exclude ascent from take-off and descent to landing as well as to avoid 
turns that were very close to either end point. Once these times were determined, 
the individual measurements were selected using a heuristic randomised process. 
The intent of this process was to avoid manual selection bias and to create measure- 
ment sets that emulate the data available for the accident flight. Measurements were 
selected recursively. 

Let t,..; denote the measurement time for the previous measurement; to is the man- 
ually selected starting time. Each measurement has a collection time labelled t;, for 
j € (I, ..., J}, where J is the total number of measurements in the communication 
log. The first measurement was selected by assigning a probability 


F 1 
p;(0) = P(0) 'ew[-5; (j - 9]. (9.1) 
J 1 j 
P(0) -Xee[-u — to) | (9.2) 


where o was chosen to be 15 min. The selected measurement was then chosen by 
taking a single multinomial draw on the probability vector p (0). This selection prefers 
measurements closer to the start time. Subsequent measurements were chosen with 
a mean time spacing of 1h. Let /(7) index the measurement chosen as the ith in the 
sequence. A probability vector for the (i + 1)th measurement was defined as 


PG + 1)7 exp [-z [b des y'] 2 jd 
0, j xl), 
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Pi+l= 5 exp | 552 (tj fi) y] : (9.4) 
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Fig. 9.13 Example measurement timings for a single flight 


Measurement (i + 1) is again selected using a single draw on a multinomial distrib- 
ution defined by p;(i + 1). The process concludes when the measurement selected 
occurs after the desired end time: this measurement is discarded. 

Figure 9.13 shows an example of the measurement times for the ten different sets 
generated for the Kuala Lumpur to Amsterdam flight on 26 February 2014. Squares 
are used to mark the initialisation time and the final time point, neither of which 
have measurements. The measurement times are marked with circles. Each row is 
a realisation of the measurement selection process. Some measurements are used 
by more than one of the sets. The number of measurements selected varies between 
eight and ten, the duration of the flight segment is approximately seven hours and 
35 min: seven one-hour spaces would lead to eight measurements in seven hours. 


9.7.2 Performance Measure 


In the object tracking literature it is common to use accuracy measures to quan- 
tify tracking performance, for example [7]. Accuracy measures quantify how well 
the estimates from the tracker match the truth. The most frequently used accuracy 
measure is root-mean-square (RMS) error, which is typically the average geometric 
distance between the true object position and the tracker estimated position. The 
requirement for MH370 is a search region, not a point estimate, so RMS is not 
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applicable. The other common accuracy measure is the Normalised Estimation Error 
Squared (NEES). This is defined as the expectation of the inner product of the esti- 
mation error with itself, normalised by the estimator covariance. For a scalar, this 
is the mean squared error divided by the filter covariance estimate. Whereas RMS 
quantifies how accurately the filter finds the centre of mass of a distribution, NEES 
quantifies how accurately the filter estimates the spread of a distribution. NEES 
inherently assumes a uni-modal distribution. Again, NEES is based on an assumed 
Gaussian system with a point estimate and covariance estimate. It is not an appro- 
priate measure for the multi-modal pdf produced by the filter in this application. 
Instead, the statistical performance of the filter output was quantified by measuring 
the highest posterior density (HPD) interval at the true aircraft location. The HPD 
interval is defined as the spatial region for which the filter output pdf is at least as high 
as the value at the true location. Figure 9.14 shows an example of this process for a 
scalar random variable x with a Gaussian mixture pdf p(x). The two components 
are equally weighted, one with mean 2 and variance 0.25 and the other with mean 5 
and unit variance. Supposing that the truth in this case was x — 6, the HPD interval 
is shaded in red and corresponds to the regions in x for which p(x) > p(6). Because 
the distribution p(x) has two modes and the value of p(6) is between the lower peak 
and the intermediate minimum, the HPD is composed of two intervals. If the truth 
had been 2.5 instead then only one region around the higher peak at 2 would be in the 
HPD and if the truth were 8 then almost all of the pdf would be in the HPD region. 

The integral of the pdf over the HPD interval corresponds to the cumulative 
probability that a random sample from the distribution is more likely than the truth 
point. If the integral is close to unity, then the HPD interval contains most of the 


p(x) 


Fig. 9.14 Highest posterior density interval of a Gaussian mixture 
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support of the pdf, that is the truth point is at a very low part of the pdf. Alternatively, 
if the integral under the HPD is close to zero then only a small portion of the event 
space is more likely than the truth point. 

Mathematically, the HPD integral is given by 


h(x" ^: p(x) = J p(x)dx, (9.5) 
x: p(x)>p(xtuth) 


where x''''^ is the true aircraft location and p(x) is the filter output pdf. In the 


discussion that follows, we abbreviate as h = h(x"; p) the random variable 
derived by transforming the random variable x" using (9.5). 

If the truth values were indeed random samples from the filter output pdfs, then it 
is relatively easy to show that the distribution of h would be uniform on the interval 
[0, 1].! If integrals tend to be clumped closer to zero then the pdfs being assessed are 
pessimistic: the tails decay too slowly and the coverage of the pdf is too broad. If the 
integrals tend to be clumped closer to unity then the truth is always in the tails and 
the pdfs being assessed are overly optimistic. For the MH370 search definition we 
prefer a conservative pdf that is a little pessimistic, in order to minimise the chance of 
excluding the true aircraft location. Provided the search zone defined can be feasibly 
measured it is better to make this region a little too large and improve the likelihood 
that the truth is contained. 

For each flight we have only ten different measurement sets so it is not feasible to 
construct a sensible estimate of p(/1). Instead we plot an estimate of its cumulative 
distribution and compare it with the line y — x, which is the cumulative distribution 
of a uniform random variable. If the h values are relatively small then the empirical 
cumulative distribution function (cdf) will rise more quickly than the reference and 
the curve will be above it. Conversely if the values are relatively large then the 
empirical cdf will rise slowly and the curve will be below the reference. 


9.7.3 Results 


Figure 9.15 shows the empirical cdf derived for each validation flight separately. This 
shows that the results within a single flight are quite correlated because the filter 
performance is dependent on geometry. For the Mumbai to Kuala Lumpur, Kuala 
Lumpur to Beijing and Beijing to Kuala Lumpur flights the h values are generally 
small but not close to zero. This indicates that the spread of the filter pdf is too high 


'To see this, let Y = p(x), i.e., the random variable obtained by applying the random value x to 
its pdf. Then the cumulative distribution function (cdf) of Y, Fy (y) = P(Y < y) is one minus the 
HPD value in (9.5). It is well-known that, assuming continuity and monotonicity of the cdf, the 
random variable obtained by passing a random value through its cdf is uniform on the interval [0, 1] 
(e.g., [13]), and if Y is uniform on [0, 1], then so is 1 — Y. The necessary assumptions are satisfied 
if the pdf p(x) contains no non-zero flat regions and no Dirac delta components. 


84 


cdf 


cdf 


cdf 


0.8 


0.6 


0.4 


02 


0.8 


0.6 


0.4 


0.2 


0.8 


0.6 


0.4 


0.2 


Kuala Lumpur — Amsterdam 26-2 


0.2 0.4 0.6 0.8 
HP integral 


Kuala Lumpur - Bejing 


0.2 0.4 0.6 0.8 
HP integral 


Kuala Lumpur — Amsterdam 7-3 


* 


0.2 0.4 0.6 0.8 
HP integral 


cdf 


cdf 


cdf 


0.8 


0.6 


0.4 


0.2 


0.8 


0.6 


0.4 


0.2 


0.8 


0.6 


0.4 


0.2 


9 Validation Experiments 


Mumbai — Kuala Lumpur 


0.2 0.4 0.6 0.8 
HP integral 


Beijing Kuala Lumpur 


* 


0.2 0.4 0.6 0.8 
HP integral 


Kuala Lumpur — Frankfurt 


X 


x 


0.2 0.4 0.6 0.8 
HP integral 


Fig.9.15 Cumulative density plots for individual validation flights. Crosses show individual data- 
set results and the solid line shows the theoretical result for independent samples 


9.7 Quantitative Analysis 85 


and that the peak of the pdf is biased. The bias occurs because the flights make small 
manoeuvres that are unobservable by the filter. For example, Fig. 9.5 shows that in the 
Kuala Lumpur to Beijing flight the aircraft made a number of heading changes that 
lasted for only a short time before the heading reverted back to its previous value. The 
minor course corrections result in a displacement in position. The filter will sample 
these paths but their dynamics are less likely than paths without a manoeuvre. For 
these flights the mode is not a reliable indicator of the true aircraft location but a 
fairly tight interval is. 

In the longer Asia to Europe flights the A values tend to spread between 0.25 and 
0.8. Again there is bias due to the repeated geometry and very large values are not 
observed because the model allows manoeuvres that are more dynamic than what 
occurred in the actual flights and this spreads the pdf. 

Figure 9.16 combines all of the trials into a single A cumulative distribution. In 
this plot the two different groups of flights are apparent: there is an initial very sharp 
rise due to the contributions of the intra-Asia flights and then a gradual climb from 
the Asia to Europe flights. 

Overall the results show that for all of the flights and measurement combina- 
tions tested the true aircraft location was inside a 85% confidence region of the 
pdf. That is, the largest h value observed was approximately 0.85. This means that 
the pdf estimates are conservative. The spread of the estimated pdf is wider than 
the spread of true values. This occurs for two reasons: firstly, the aircraft dynamic 
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Fig. 9.16 Cumulative density plot combined over all validation flights. Crosses show individual 
data-set results and the solid line shows the theoretical result for independent samples 
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model provides more flexibility than is typically used; for example, in most com- 
mercial flights, smaller turns are more likely than turns of 90? or more. Secondly, the 
assumed measurement variances were deliberately inflated to be pessimistic, as dis- 
cussed in Sect. 5.3. Given that the accident flight was not a typical commercial flight, 
the dynamic model should not be exactly matched to typical commercial flights. A 
somewhat conservative pdf in this case is desirable so long as the pdf does not spread 
over an area that is unreasonably large to search. 
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Chapter 10 
Application to the MH370 Accident 


The previous chapters have constructed a Bayesian method for estimating 
commercial aircraft trajectories using models of the information contained in satel- 
lite communications messages and of the aircraft dynamics. This chapter applies the 
estimator to the accident flight. We show the pdf at the final message time 00:19 and 
perform several tests that provide confidence in the robustness to the model parame- 
ters and measurement characterisation. Probable end of flight dynamics are used to 
convert the pdf in the air at 00:19 into an ocean search zone. As a further measure of 
robustness, we show that the model can be initialised at the end of the initial ascent, 
as was done in the validation flights, and the output pdf remains in the same broad 
geographic region, even without using the subsequent primary and secondary radar 
data. 

Table 10.1 summarises the available measurement data and the BTO measurement 
standard deviation, when applicable, since this depends on the message type. The two 
anomalous R1200 messages and the R600 message occurred during transient phases 
of operation for the SATCOM equipment so the BFO values reported for these times 
cannot be used. The BTO measurement variance is also inflated to account for the 
uncertainty in the correction term. 


10.1 The Filter Applied to the Accident Flight 


This section summarises the results of applying the variable rate particle filter to the 
data available for the accident flight. We consider two different treatments of the 
measurements. Figure 10.1 shows the Google Earth pdf representation for the case 
where only BTO measurements were used. As with the validation flights, the colour 
coded paths show the most likely routes to each latitude. Clearly this distribution 
is multimodal. Based on BTO data, it can be deduced that the aircraft must have 
manoeuvred in some way between 18:28 and 19:41. However the BTO measurements 
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Table 10.1 SATCOM messages used for MH370 analysis 


Time (UTC) Measurement type | BTO BTO std. dev. (us) | BFO 
18:25:34 Anomalous R1200 | Y 43 N* 
18:28:05, 18:28:14 | R1200 Y 29 Y 
18:39:55 C-channel N - Y 
19:41:02, 20:41:04, | R1200 Y 29 Y 
21:41:26, 22:41:21 

23:15:02 C-channel N - Y 
00:10:59 R1200 Y 29 Y 
00:19:29 R600 Y 63 N* 
00:19:37 Anomalous R1200 | Y 43 N* 


Measurements marked with an asterisk are available but cannot be used as discussed in the text 


Ggoale earth 
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Fig. 10.1 Set of paths from 18:02 to 00:19 using only BTO measurement weighting (i.e., not using 
any BFO measurements) 


are consistent with two options: heading broadly North or broadly South. Many of 
the paths turn very early, almost immediately after initialisation. This is because the 
BTO arc at 18:25 is not consistent with the radial speed prescribed by the initial 
speed and direction. The 18:22 radar point, while potentially inaccurate, implies that 
no turn had occurred at that stage and that the speed may have changed. Due to 
the questions surrounding its accuracy, the filter does not use this radar point, so it 
postulates early direction changes as well as early speed changes to match the 18:25 
BTO. A feature of most of the paths is that they do not choose to make subsequent 
turns even though there is no BFO data to influence angle likelihood. The BTO data 
itself from 19:41 onwards is consistent with straight and level flight. 
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Fig. 10.2. Set of paths from 18:02 to 00:19 using BTO and BFO measurement weighting. Purple 
box shows the November 2015 indicative search region 


Figure shows how this path set is pruned by including BFO data. The crucial 
BFO measurements are at 18:28 and 18:39. The first of these is consistent with the 
same heading as at initialisation, so it implies that the aircraft has not yet turned. 
In contrast, as discussed in Sect. 5.4, the 18:39 BFO measurement indicates that the 
aircraft is moving roughly South at this time. Between the two, these measurements 
restrict the time of the turn to a window between 18:28 and 18:39. This is reflected 
in the set of surviving paths. Paths that went very far South had to turn earlier and 
have been rejected as have all of the paths heading to the North towards Asia. The 
plot also shows the indicative search region as at November 2015, displayed as a 
purple box. 

The main effect of including BFO data is to resolve the ambiguity about the 
manoeuvre after 18:28. The BFO data does not significantly change the shape of the 
Southern mode of the pdf. To show this more clearly Fig. plots the two pdfs, 
without and with BFO data, as one dimensional curves parameterised by latitude. 
The latitude is a nonlinear function of the position around the BTO arc so these pdfs 
are distorted slightly, but they clearly show the effect of the BFO data: it selects the 
mode from a multimodal pdf but it does not significantly change the mode shape. 
Because of this, including only the C-channel BFO measurement has the same effect 
and changing the assumed BFO noise standard deviation has no significant effect 
either unless it could be artificially reduced to less than 1 Hz. We conducted several 
experiments that varied the assumed BFO noise and which BFO values were included. 
The results are consistent with the curves in Fig. 1C 
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Fig. 10.3 Comparison of the filter pdf without (top) and with (bottom) BFO measurements 


10.2 Manoeuvre Statistics 


An interesting question is whether the posterior pdf output from the filter provides 
any strong estimate of the number and type of manoeuvres after the loss of primary 
radar coverage. To address this we compare the probability mass of the number of 
turns for the case of no measurements, i.e., the prior probability from predictions 
used to create the enormous disc pdf in Fig. 7.4, and the posterior probability after 
applying the measurements. These are shown in Figs. 7.3 and 10.4 respectively. The 
posterior distribution does not count turns made between 00:11 and 00:19 because 
the BFO data at 00:19 is not used (as it is believed to be unreliable); this permits 
the filter to make superfluous manoeuvres. The histograms were created by counting 
the number of manoeuvres performed by each particle and then scaling their con- 
tribution to the histogram by the particle weight. In the prior probability there is a 
significant probability of making a very large number of manoeuvres. The particles 
that contribute to this part of the histogram are those that have selected very short 
mean manoeuvre times T. 

The histograms show that the posterior number of turns has reduced dramatically, 
as has the number of speed changes. The likelihood of selecting a sequence of random 
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Fig. 10.4 Histogram of posterior number of turns and speed changes 


speeds that average out to a suitable overall speed is relatively high. Since there 
is usually at least an hour between BTO arcs, the filter has time to make several 
speed changes. Only the average of these speeds influences the position and hence 
the BTO. As discussed several times, the BFO is sufficiently uninformative that the 
instantaneous velocity at the measurement time is only loosely constrained. However, 
itis unlikely that a sequence of random turns will “cancel each other out,” and the filter 
finds very few paths that string together multiple turns when a straight path would 
suffice. The low support for turns leads to selection of longer mean manoeuvre times, 
which in turn suppresses superfluous speed changes. 

The turn histogram indicates that around half of the paths made more than a single 
turn. This would appear to be of interest but is in fact misleading. Figure 10.5 shows 
the posterior distribution of angle as a function of time. The greyscale image shows 
the probability density at each time slice; the darkest points are most likely. The solid 
line shows the mean and one-sigma error bars indicate the standard deviation. Note 
that the horizontal scale is measurement index, not time, but time values are used 
to label the measurements to give perspective. The angle plot shows that there are 
a number of particles that make a very small turn before 18:28 and then follow this 
with a second turn between 18:28 and 18:39. The result is very similar to making a 
single turn through the same total angle. Similarly, there are a number of particles that 
make a turn before 18:39 and then a second very small turn immediately afterwards. 
Again, this is approximately equivalent to a single turn through the same total angle. 
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Fig. 10.5 Posterior control 320 
angle distribution 
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The apparent large number of double turns is really just a single turn being broken 
into two segments. There appears to be very few genuine turns later in the flight. 

Figure 10.6 shows the marginal distribution of the Mach number between 18:02 
and 00:19. The prior is uniform in Mach between 0.73 and 0.84 but the posterior 
strongly prefers the higher speed part of this range. This is approximately consistent 
with the most likely speeds for the aircraft over long durations. 
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Fig. 10.6 Posterior speed distribution 
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10.3 Residual Measurement Errors 


Another measure of the quality of an estimate is the statistical behaviour of the 
measurement residual errors. If the estimates were exactly the truth then the residual 
measurement errors would be the actual measurement noise. So these residuals should 
be uncorrelated, zero mean, and spread consistently with the measurement noise 
variance. Figure 10.7 shows the marginal distributions of the residual measurement 
errors. The filter predicts the BTO and BFO at the measurement times based on 
the estimated position and velocity and the residual is the difference between these 
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Fig. 10.7 Residual measurement error 
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predictions and the measured values. Each sample is weighted according to its overall 
likelihood and the combination gives the marginals in Fig. 10.7. The horizontal axis 
is again the measurement index and the vertical axis is residual error, dark points 
have high probability, and the error bars show the one-sigma limits. The figure also 
overlays the one-sigma lines for the assumed measurement noise. For BTO this varies 
with the message type, for BFO it is a fixed +7 Hz. There is no BTO measurement 
for the C channel communications that occurred at 18:39 and 23:14. The BFO has 
been discarded at 18:25 and 00:19 since these readings are thought to be unreliable. 
The plots only show residuals when a valid measurement is available. 

For the BTO residual there are two larger residuals on the final two measurements. 
These measurements are actually very close together in time: the first is an R600 
message at 00:19:29 and the second is an R1200 message at 00:19:37. The reason 
why both of these measurements show large residuals is that they are not consistent 
with each other: the residuals have opposite signs, reflecting that one measurement 
is longer in range than the prediction and the other is shorter. In the absence of a 
reason to prefer one over the other, we use both measurements and let the filter find 
paths that are the best statistical fit. The BFO residuals are statistically consistent 
with the empirical error model. 


10.4 Posterior Distribution of Manoeuvre Time Constant 


The inference procedure that has been developed samples the manoeuvre time con- 
stant, which is a static parameter. In the case of conventional particle resampling, this 
is not advisable, as subsequent resampling steps generally leave few distinct values 
of the static parameter. The recursive resampling process described in Chap. 8 is less 
prone to this difficulty, as many independent sample paths are retained. 

The sample support of the posterior distribution of the manoeuvre time constant 
is demonstrated in Fig. 10.8. The diagram clearly shows that the distribution is non- 
degenerate. The figure compares the posterior distribution from the filter, which is 
shown as bars, with the Jeffreys prior, which is shown as a red line. The Jeffreys 
prior prefers lower mean manoeuvre times but the data does not support these. The 
posterior is dominated by longer mean manoeuvre times, with 97 96 of the distribution 
having a mean time between manoeuvres of more than one hour, and 83 % having a 
mean time more than two hours. 


10.5 End of Flight 


The output of the particle filter is an estimate of the pdf of the aircraft state at 00:19. 
The aircraft was still in the air at this time and a model is required to describe the 
distribution of how it may have descended. This has been primarily the responsibility 
of the ATSB and the other members of the search team. A discussion of the different 
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methods used to model the potential motion is presented in [5]. The model for aircraft 
motion after 00:19 leads to a prioritisation of the search along and around the final 
BTO arc. 

Quantitatively, the model of descent defines a transition distribution p(Xtinai|Xx ), 
which describes the probability density of the final state given the state at time 00:19. 
In effect this acts as a kernel to spread the distribution at 00:19 via the expression 


P P 
p (XtinallZx) = J P(Xtinall XK) X wiS(Xy — xk)dxg = S wk p(Xtinail Xe) 
=1 =1 

” ” (10.1) 
The analysis in [5] leads to a probable scenario where the aircraft ran out of fuel 
at some time between 00:11 and 00:19. The final satellite communications message 
could be due to the modem rebooting under auxiliary power. Under this hypothesis, 
the aircraft engines were already unpowered at 00:19. The spread of the kernel 
function is then determined by the distance over which the aircraft could have moved, 
which depends on whether or not the aircraft was under human control during this 
period. 

Flight simulator studies of uncontrolled descents have shown a high likelihood 
of the aircraft reaching zero altitude within 15 nm of the beginning of descent [5]. 
However, the beginning of descent is not known. It is possible for the aircraft to have 
travelled farther, especially if a human was controlling the aircraft. As an indicative 
kernel, and following advice from the ATSB, a uniform disc of radius 15nm with 
a Gaussian drop off with standard deviation 30nm beyond this was chosen; this 
represents the accident investigators' assessment of the likely scenarios. Figure 10.9 
shows a radial slice through this kernel function. The kernel was convolved with the 
particle locations at 00:19 to generate the heat map shown in Fig. 10.10. The black 
region overlaid on the heat-map shows an indicative bounding box on the geographic 
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Fig. 10.9 Descent kernel 


kernel 
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search region. The figure also shows the 99, 95 and 90 % confidence regions, i.e., the 
smallest geographic area containing each percentage of the pdf. 


10.6 Earlier Initialisation 


As a further level of validation on robustness, a separate run of the filter was performed 
initialising using the final ACARS report. This run is not useful for defining a search 
region because it ignores important data from the primary and secondary radar. 
However, the result does illustrate the ability of the method to address flight paths with 
several large turns. Figure 10.11 shows the pdf derived from the 17:01 initialisation. 
The pdf is overlaid with the one that is presented earlier in this chapter that uses all 
of the available information. Two versions of the 17:01 initialisation are presented. 
The first, which ignores the primary radar data, is shown in red in the upper diagram. 
The second, shown in the lower diagram, re-weights the paths by treating the 18:02 
radar point as a measurement with 1? standard deviation in latitude and longitude. 
This weighting is rather loose but the samples are limited and none of them are 
sufficiently close to the radar point to permit a tighter match. The filter has given a 
similar output even though the initialisation point was heading North and there are 
two very significant manoeuvres reported by the radar in the first hour. This provides 
further confidence in the ability of the filter to characterise aircraft motion even in 
complicated situations. The sequence of turns performed by MH370 between 17:01 
and 18:02 is much more dynamic than any experienced in the validation flights. 
When the radar data was ignored the filter found a cluster of paths that travel a 
significant distance to the North before turning prior to 18:39. These paths conclude 
atthe Northern edge of the search region because they travel too far to the North early 
on, in a manner which is inconsistent with the primary radar data. In contrast, the 
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SN Google earth 


17:01 initialisation weighted by proximity to 18:02 radar (red) compared with the 18:02 radar based 
initialisation (white) 


Fig. 10.11 Probability density from early initialisation 


paths that are preferred after re-weighting using the primary radar data travel South 
of the radar point and arrive at the Southern edge of the search region. One would 
expect that, given further computation, samples would be found between these two 
that would match the 18:02 radar point better and agree better with the search region. 
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10.7 Cost Index 


Another autopilot mode for control of air speed is Cost Index, as discussed in 
Sect. 6.2.1. To recap, under this mode the autopilot dynamically selects a speed in 
order to optimise a criterion which trades fuel consumption and travel time accord- 
ing to the Cost Index value entered by the pilot [34]. An additional experiment was 
performed in which the aircraft was permitted to enter into the Cost Index mode 
at a random time, remaining in that mode from that time forward. The Cost Index 
value was randomly sampled between 0 and 100, and lookup tables provided by 
Boeing were used to determine a speed based on the sampled value of altitude, and 
an estimate of the aircraft weight, which varies over time as fuel is expended. 

The result of this experiment is shown in Fig. 10.12. The white curve shows the 
pdf calculated through the model using Mach number, as shown in Fig. 10.2, while 
the red curve shows the pdf calculated using the Cost Index model. The figure shows 
that the use of Cost Index does not significantly change the result; the distribution is 
slightly more compact than in the case of the constant Mach number model. Since 
either mode is possible, the broader result in Fig. 10.2 is preferred for determining 
the search region. 


Google earth 
€ 


Fig. 10.12 Set of paths obtained using Cost Index mode for controlling air speed (red), compared 
to result using the Mach number model (white) as shown in Fig. 10.2 
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10.8 Other Variations 


The results presented in this chapter are a summary of the analysis we have undertaken 
on the accident flight data. They are far from an exhaustive set of the variations that 
have been examined. Many other changes to the filter parameters and measurement 
model were explored but none made a significant impact on the search region. Exam- 
ple changes to the filter include: modifying the Jeffreys prior; using fixed average 
manoeuvre times; changing the incremental time step; resampling the Jefferys prior 
to emulate a change in aircraft behaviour; moving the initialisation point; modifying 
the initialisation variances; forcing a prescribed angle after 18:39; and quantising 
heading to integer degrees. Example changes to the measurement model include: 
vetoing selected measurements; using only C-channel BFO measurements; increas- 
ing and decreasing the BFO assumed noise variance; replacing the BFO with a 
"true" BFO using known location for the validation flights; including or excluding 
particular BTO measurements and R600 messages; and cross validation where each 
measurement is removed in turn to investigate output robustness. 

In our earliest attempts at producing a pdf a much simpler aircraft model was 
used. For example, under the assumption of a single turn, the parameter vector is 
relatively short and more conventional estimation methods can be used. However, 
such a model cannot be applied to any of the validation flights and runs the risk of 
constraining the estimated path with additional assumptions. The fact that the very 
general model here chooses to make few turns provides information that would not 
be available from a single turn geometric model. 

The factors that do make a significant difference to the output pdf are the assumed 
spread of Mach number and the end of flight model. The assumed Mach number 
range covers the speeds feasible to achieve the required flight endurance time. The 
lower end of this speed range results in the Northern part of the pdf and the higher 
end of the speed range results in the Southern part. Restricting the speed to only 
Mach numbers above 0.8, for example, would contract the pdf to the South. The 
consequence is that using a smaller speed range within the bounds already modelled 
leads to a subset of the search zone. If a different end of flight model is assumed 
the general consequence is to spread the search zone over a larger area. Simulations 
have predicted that the maximum distance that the aircraft could have glided under 
human control is approximately 100 nm after 00:19 [5]. The search zone that this 
scenario would imply is very much larger. 
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Chapter 11 
Ongoing Refinement 


The previous chapters of this book have built up models for the limited available 
measurement data and for aircraft motion, and then used these models to produce 
a pdf of the final aircraft latitude and longitude. In principle, one could use this 
pdf to direct a search and then the very act of searching would provide further 
measurement data. This chapter describes how the Bayesian method can be used to 
adapt the position pdf to account for data collected after the accident flight. Two 
data sources are discussed: the sonar imagery data collected in the search and the 
discovery of the flaperon on Reunion Island. 


11.4 Updating the Distribution Using Search Results 


The measurements collected as part of the search can be treated in the same mathe- 
matical framework as the communications data. This was used, for example, in the 
search for AF447 [40]. In this case the aircraft is no longer moving so the prediction 
stage becomes degenerate and the predicted pdf is the same as the previous posterior 
pdf. The result of the search can be summarised by the probability that the cumulative 
search effort would have detected the aircraft at any particular location, Pp (x). The 
posterior pdf given the search effort (denoted S) then becomes: 


PXtinallS, ZK) X [1 — Pp(Xtina)] P (XtinallZi) , (11.1) 


where the constant of proportionality is determined to ensure that the posterior inte- 
grates to unity. If a particular area A is searched with a constant probability of 
detection Pp, the probability of finding the aircraft is 


P (find during search of A) = Pp I P (Xtinall Zik) dXtinal- (11.2) 
A 
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In the analysis used for prioritising the search for flight AF447, the probability of 
detection for areas searched using side-scan sonar was modelled as 0.9 [40]. This 
type of analysis could be used, for example, to prioritise whether to revisit a more 
likely area or to search a less likely area that has not yet been searched. 

Based on the quality assurance process which has been implemented in the MH370 
search, which includes revisiting items that have been assessed as potential debris, 
it is considered highly unlikely that the search would fail to detect the aircraft if 
the correct location is searched. Due to sensor drop-out and terrain masking, there 
will inevitably be small pockets which are not covered in a first pass of the search, 
and (11.1) can be used to determine the priority of returning to ensure that these are 
examined. 


11.2 Reunion Island Debris 


The apparent lack of debris from the aircraft was a mystery which was resolved in 
part when a flaperon was discovered washed up on Reunion Island on 29 July 2015, 
508 days after the accident, and later confirmed to be from MH370. Two questions 
arise from this find: 


1. What information does the discovery of the flaperon on Reunion Island provide 
about the final location of the aircraft? 

2. What information does the lack of any other debris to date provide about the final 
location of the aircraft? 


Debris information was used to inform the search for AF447 [40], and was pro- 
posed for MH370 in the white paper [19]. The most directly relevant source of 
information on the likely drift of the flaperon is from the Global Drifter Program 
[32]. The program provides 30 years of data on buoys referred to as drifters, which 
are regularly deployed in oceans worldwide. The primary goal of the program is to 
measure ocean currents rather than surface effects. To ensure that the drifter motion 
is dominated by the ocean current, the drifters are deployed with drogues, which are 
sea anchors that sit around fifteen metres beneath the ocean surface. If the drogue 
is detached then there is an observable change in the motion of the drifter [28]. A 
drifter without a drogue is referred to as an undrogued drifter and is more buoyant 
than the drifters with an attached drogue. The buoyancy of the flaperon has not been 
characterised at this time but it is expected that the undrogued drifters better describe 
its likely motion because the drogues by design cause drifters to move according to 
deeper currents. 

In this section, we derive an updated pdf of the final location of the aircraft based 
on an analysis of data from the Global Drifter Program. It should be noted that the 
characteristics of the drifters are not accurately matched to the flaperon, thus there 
is some uncertainty in the applicability of the results, but it is believed to be the 
most relevant data source available. The authors gratefully acknowledge Dr. David 
Griffin (CSIRO, Australia) for providing specialist analysis and advice in support of 
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this effort, as well as the results illustrated in Fig. 11.1. Further information may be 
found at [17]. 

The first question is answered by modifying the pdf provided in the previous 
chapter. This is achieved by developing a likelihood function of the source location 
of the flaperon using the drifter data. The likelihood is used to calculate the posterior 
pdf of the final aircraft location given both the information provided by the Inmarsat 
satellite data and the discovery of the debris. 

The second question relates to the potential information provided by the lack of 
other debris. We consider a statistical framework for quantifying this information 
but the framework requires parameters that cannot be reliably determined. For this 
reason we limit ourselves to qualitative conclusions about the lack of other debris. 


11.2.1 Update of Final Location Distribution 


We start with the form which results from the particle filter: (sampling to incorporate 
the effect of the descent kernel of Sect. 10.5) 


P 
P Gcinail Zi) © >) w^ (Xrinai — x^) , (11.3) 


p=1 


For a single item of debris, we seek a model which permits us to update (11.3) using 
the knowledge that the item arrived at a given location and time. For this, we require 


a transition distribution for how the debris would move over 508 days. Denoting the 
location after 508 days as y, this transition distribution is p(y|Xtinai). Thus, assuming 
a single item of debris, the kinematic distribution updated by the knowledge of the 
debris find at location y = y is: 


PtinallZk. Y = y) X P(Y = ylxtina) P Xtinall Zr) 


P 
~ X w”p (y = yix”) 8 (Xma —x"), — (14) 
p-l 


Equation (11.4) can be interpreted as re- weighting each particle x? by the likelihood 
that an item of debris at that location would end up at location y — y. 


A field of debris can be modelled as a Poisson point process (e.g., [42]). Suppose 
that the expected number of debris items is å, that an item that is washed up on the 
shore is identified with probability P;, and that no items are found while they remain 
in the ocean. Then the likelihood of the first item of debris being found at y — y 
given that the end-of-flight location was Xfinal is: 


I(debris|xina) = exp [-m pirixwaDdy 
C 


x exp (—Ap(y = ylXtina)) Ap(y = ylXtina) (11.5) 
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where the integral is over the coastal region C, i.e., the region where debris would be 
likely to have been found. This expression could be used to update the model as: 


P Xtinai|Z.« , debris) ox I(debris[xrnai) p (Xtinai| Z.« ) 
P 
& S w”l (debris|x”) à (xinai — x”) , (11.6) 


p-l 


Similar to the case of a single debris item, this model re-weights each particle x?” with 
the likelihood that a debris field starting from that location would not yield any identi- 
fied coastal debris within the time other than a single item at y — y. The model could 
be easily refined to incorporate thinning, i.e., the likelihood that over time, items of 
debris might sink, as well as the temporal aspect that the debris did not arrive until 
some 508 days later, by splitting the transition distribution p(y|x) into a number of 
intermediate, smaller duration transition steps. However, the required model para- 
meters, are very difficult to quantify with any degree of confidence. Specifically, the 
expected number of debris items, and the probability that an item washed up would 
be identified, are assumed to be known. Since these parameters cannot be adequately 
specified, we rely on a qualitative approach for incorporating information provided 
by the lack of debris, and use the quantitative approach of (11.4) for incorporating 
the information on the debris item that was found. 


11.2.2 Data from Global Drifter Program 


The data available for undrogued drifters are summarised in Fig. 11.1. Each diagram 
shows drifter trajectories which pass through a region of interest in months between 
February and April (i.e., the time of year of the accident). The middle figure centres 
the region of interest on the search zone, while the top figure examines the region 
to the North, and the bottom figure shows the region to the South. Coloured dots 
are drawn along each drifter trajectory every 100 days to illustrate time along the 
trajectory. Trajectory segments coloured in blue indicate water colder than 18?C, 
the temperature above which barnacle nauplii settlement and growth is accelerated; 
evidence of accelerated growth was present on the flaperon [17]. 
The following observations can be made from the figures: 


1. In the top diagram (the region to the North of the search zone), many trajectories 
head West at quite a fast rate, reaching Madagascar and Tanzania in 300 days or 
less. The absence of debris being identified in the Western Indian Ocean many 
months earlier than July 2015 would tend to indicate this region as being less 
likely. 

2. Asignificant proportion of the trajectories in the middle diagram (the search zone) 
arrive in the general vicinity of Reunion Island at around 500 days. 
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Tracks from Feb-April entry into blue rectangle to 508days (max) later. Blue tracks are in water colder than 18C. 
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Fig. 11.1 Trajectories from global drifter program. Top diagram shows trajectories which pass 
through the box with latitudes from 27°S to 33°S and longitudes from 96°E to 106?E. Middle 
diagram shows trajectories which pass through the box from 33 °S to 39 °S and longitudes from 87 °E 
to 96°E. Bottom diagram shows trajectories which pass through the box from 39°S to 45°S and 
longitudes from 77 °E to 87°E. Coloured dots show time, marking every 100 days. Blue trajectories 
are in water colder than 18°C, the temperature above which barnacle nauplii settlement and growth 


is accelerated. Figure courtesy of Dr. David Griffin, CSIRO; see also [17] 


106 11 Ongoing Refinement 


3. The majority of trajectories in the bottom diagram (the region to the South of the 
search zone) head to the East, towards Australia and New Zealand. Very few head 
toward the vicinity of Reunion Island. 


Thus the discovery of the single flaperon on Reunion Island after 508 days would 
seem to be more consistent with the existing search region, and less consistent with 
regions significantly to the North or South. 


11.2.3 Posterior Distribution Using Debris Data 


Figure 1 1.1 represents the total undrogued drifter data available passing through each 
box at the relevant time of the year over the 30 year history of the Global Drifter 
Program. It is insufficient to construct a transition distribution, hence additional 
processing was performed to enhance the data set. Specifically, pairs of trajectories 
that passed close together (possibly in different years) were “joined,” i.e., used to 
create additional synthetic trajectories comprised of the head of one with the tail of 
the other. A synthetic data set was generated by joining four trajectory segments. In 
each joining process, each previous trajectory creates around 30 new trajectories. 
The joined trajectories provide a sampled representation of the distribution: 


P 
1 
Pariter X, y) © 5 9 ,8(x — xP) ly — y^) (11.7) 


p-l 


where x is the location of the drifter on 8 March 2014, and y is the location on 
29 July 2015. As in (11.4), the distribution of interest is p(y|Xtina) evaluated for y 
at Reunion Island, or, pragmatically, integrating over a small region R surrounding 
Reunion Island: 


I (Xtinai) zj P(y|Xtinal)dy (11.8) 
R 


The conditional distribution p(y|x) requires normalisation by the prior distribution 
Parifter(X); this represents the overall mean density of drifters in the ocean, which is 
neither spatially uniform nor time invariant. Applying a small kernel K (x) around 
each sample, the likelihood is evaluated as: 


P 
m Dárifter(Xtinai, Y) dy a P ym K (Xtina — X") + € 


Parifter (Xtinal) t x. a K (Xtina — XP) + € 


l(Xtina) = (11.9) 


where {x!,...,x?} is the set of locations for trajectories which pass though the 
region R in the right time window, including those synthetically augmented using 
the aforementioned joining process. {x!, ... , X") is the overall set of drifter locations 
for the time interval from February to April. The ¢ terms are added for regularisation, 
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Fig. 11.2 Distributions obtained using undrogued drifter data from Global Drifter Program. Top- 
left diagram shows the overall density of drifters in the months of interest (over all years). Top-right 
shows distribution of synthetically augmented trajectories which pass though the Reunion Island 
region after 508 days. Lower-left diagram shows the quotient of these two distributions, which is 
the likelihood p(y Xtinal). Lower-right diagram shows this same likelihood, zoomed in to the region 
of interest for the posterior distribution. An indicative search area is marked as a black rectangle. 
Data courtesy of Dr. David Griffin, CSIRO; see also [17] 


i.e., preventing very large values in the conditional resulting from the finite sample 
support of the denominator (we set € = 107^). 

The result of this analysis is shown in Fig. 11.2, using a fixed kernel of standard 
deviation 1? in latitude and longitude, and evaluating the density estimates using 
the KDE toolbox [22]. The top left diagram shows the overall spatial distribution of 
drifters; the indicative search area is again shown in black. The overall distribution 
clearly has significant spatial variation, particularly in the search area. The density of 
drifters that pass though the Reunion Island region after the right duration is shown 
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Fig. 11.3 Posterior distribution based on Inmarsat data (top), and incorporating debris discovery 
at Reunion Island (bottom) 


in the top right. It exhibits spatial variations very similar to the overall variations 
in the top left. The likelihood function p(y|Xtinai) is the ratio of these two densities 
and is shown in the lower left. Note that large values can occur in areas of limited 
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sample support due to noise. The variation in the distribution of drifters that arrived at 
Reunion Island is almost completely explained by the variation in the overall starting 
distribution, therefore the likelihood shows much less spatial variation. The lower 
right diagram shows a zoomed in region around the search area. 

Figure 11.3 shows the posterior distribution based on the Inmarsat data from 
Fig. 10.10, alongside the updated distribution, which incorporates the information 
from the debris item located at Reunion Island. The updated distribution is shifted 
very slightly to the North, but the effect is negligible. This result is not unexpected 
given the long duration between the accident and the debris discovery. 
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Chapter 12 
Conclusions 


In this book we have described our Bayesian approach to defining the MH370 search 
zone. The three ingredients required for the Bayesian approach are 


e a prior (defined by the Malaysian military radar) 

e alikelihood function describing the relationship between BTO and BFO measure- 
ments and the aircraft state vector and 

e a model of the aircraft dynamics. 


All three are described in detail. To increase confidence in the models and process 
developed, validation has been performed using previous flights of the accident air- 
craft, and data from other flights in the air at the same time as the accident flight. 
In all cases the true location aligns with that predicted by the Bayesian analysis. 
All validation trajectories have significant numbers of speed, heading and altitude 
changes, which are successfully captured by the model. This is in stark contrast with 
the accident flight which results in a prediction of no significant manoeuvre after 
the Southerly turn near the Northern tip of Sumatra before continuing in a Southerly 
direction until it ran out of fuel in the Southern Indian Ocean, West of Australia. 

The search zone is dependent on the surface area covered by expected descent 
scenarios from the time of the final satellite log-on attempt at 00:19. This has been 
defined by expert accident investigators at the ATSB. If the actual descent scenario 
was inconsistent with the distribution of possibilities considered then the search zone 
may need to increase in area. 
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